Commit e1627801 authored by Chris Cantwell's avatar Chris Cantwell
Browse files

Merge branch 'master' into feature/dealiasing

parents 98cf6a0d 94ff5684
......@@ -8,10 +8,10 @@
</files>
<metrics>
<metric type="L2" id="1">
<value tolerance="1e-12">0.000418001</value>
<value tolerance="1e-12">0.000416575</value>
</metric>
<metric type="Linf" id="2">
<value tolerance="1e-12">0.000874515</value>
<value tolerance="1e-12">0.000871589</value>
</metric>
</metrics>
</test>
......
......@@ -8,10 +8,10 @@
</files>
<metrics>
<metric type="L2" id="1">
<value tolerance="1e-8">0.000418001</value>
<value tolerance="1e-8">0.000416575</value>
</metric>
<metric type="Linf" id="2">
<value tolerance="1e-8">0.000874515</value>
<value tolerance="1e-8">0.000871589</value>
</metric>
</metrics>
</test>
......
......@@ -9,10 +9,10 @@
</files>
<metrics>
<metric type="L2" id="1">
<value tolerance="1e-12">0.000418001</value>
<value tolerance="1e-12">0.000416575</value>
</metric>
<metric type="Linf" id="2">
<value tolerance="1e-12">0.000874515</value>
<value tolerance="1e-12">0.000871589</value>
</metric>
</metrics>
</test>
......
......@@ -9,10 +9,10 @@
</files>
<metrics>
<metric type="L2" id="1">
<value tolerance="1e-12">0.000418001</value>
<value tolerance="1e-12">0.000416575</value>
</metric>
<metric type="Linf" id="2">
<value tolerance="1e-12">0.000874515</value>
<value tolerance="1e-12">0.000871589</value>
</metric>
</metrics>
</test>
......
......@@ -8,10 +8,10 @@
</files>
<metrics>
<metric type="L2" id="1">
<value tolerance="1e-12">0.00042687</value>
<value tolerance="1e-12">0.000417674</value>
</metric>
<metric type="Linf" id="2">
<value tolerance="1e-12">0.00170172</value>
<value tolerance="1e-12">0.00168661</value>
</metric>
</metrics>
</test>
......@@ -9,10 +9,10 @@
</files>
<metrics>
<metric type="L2" id="1">
<value tolerance="1e-12">0.00042687</value>
<value tolerance="1e-12">0.000417674</value>
</metric>
<metric type="Linf" id="2">
<value tolerance="1e-12">0.00170173</value>
<value tolerance="1e-06">0.00168661</value>
</metric>
</metrics>
</test>
......@@ -392,10 +392,10 @@
<N VAR="u" VALUE="-sin(PI/2*x)*cos(PI/2*y)*sin(PI/2*z)*PI/2" />
</REGION>
<REGION REF="2">
<R VAR="u" VALUE="cos(PI/2*x)*sin(PI/2*y)*sin(PI/2*z)*PI/2+sin(PI/2*x)*sin(PI/2*y)*sin(PI/2*z)" PRIMCOEFF="1" />
<R VAR="u" VALUE="cos(PI/2*x)*sin(PI/2*y)*sin(PI/2*z)*PI/2+2*sin(PI/2*x)*sin(PI/2*y)*sin(PI/2*z)" PRIMCOEFF="2" />
</REGION>
<REGION REF="3">
<R VAR="u" VALUE="sin(PI/2*x)*cos(PI/2*y)*sin(PI/2*z)*PI/2+sin(PI/2*x)*sin(PI/2*y)*sin(PI/2*z)" PRIMCOEFF="1" />
<R VAR="u" VALUE="sin(PI/2*x)*cos(PI/2*y)*sin(PI/2*z)*PI/2+3*sin(PI/2*x)*sin(PI/2*y)*sin(PI/2*z)" PRIMCOEFF="3" />
</REGION>
<REGION REF="4">
<D VAR="u" VALUE="sin(PI/2*x)*sin(PI/2*y)*sin(PI/2*z)" />
......
......@@ -31,6 +31,7 @@ int main(int argc, char *argv[])
fprintf(stderr,"\t Chebyshev = 10\n");
fprintf(stderr,"\t Monomial = 11\n");
fprintf(stderr,"\t FourierSingleMode = 12\n");
fprintf(stderr,"\t Gauss Lagrange = 15\n");
fprintf(stderr,"Note type = 1,2,4,5 are for higher dimensional basis\n");
......@@ -81,7 +82,7 @@ int main(int argc, char *argv[])
}
else
{
Qtype = LibUtilities::eGaussLobattoLegendre;
Qtype = LibUtilities::eGaussLobattoLegendre;
}
......
......@@ -106,8 +106,8 @@ namespace Nektar
CanGetRawPtr<NekMatrix<LhsDataType, MatrixType> >
>::type* p=0)
{
int m = lhs.GetRows();
int n = lhs.GetColumns();
int m = lhs.GetRows();
int n = lhs.GetColumns();
int kl = lhs.GetNumberOfSubDiagonals();
int ku = lhs.GetNumberOfSuperDiagonals();
double alpha = lhs.Scale();
......
......@@ -469,7 +469,7 @@ namespace Nektar
{
return NekMatrix<DataType, StandardMatrixTag>(rhs, eWrapper);
}
template<typename DataType>
boost::shared_ptr<NekMatrix<DataType, StandardMatrixTag> > NekMatrix<DataType, StandardMatrixTag>::CreateWrapper(const boost::shared_ptr<NekMatrix<DataType, StandardMatrixTag> >& rhs)
{
......
......@@ -1036,7 +1036,7 @@ namespace Nektar
Array<OneD,int> sign;
StdRegions::VarCoeffMap varcoeffs;
varcoeffs[StdRegions::eVarCoeffPrimative] = primCoeffs;
varcoeffs[StdRegions::eVarCoeffMass] = primCoeffs;
LocalRegions::MatrixKey mkey(StdRegions::eMass,StdRegions::eSegment, *edgeExp, StdRegions::NullConstFactorMap, varcoeffs);
DNekScalMat &edgemat = *edgeExp->GetLocMatrix(mkey);
......@@ -1158,7 +1158,7 @@ namespace Nektar
Array<OneD,int> sign;
StdRegions::VarCoeffMap varcoeffs;
varcoeffs[StdRegions::eVarCoeffPrimative] = primCoeffs;
varcoeffs[StdRegions::eVarCoeffMass] = primCoeffs;
LocalRegions::MatrixKey mkey(StdRegions::eMass,StdRegions::eSegment, *edgeExp, StdRegions::NullConstFactorMap, varcoeffs);
DNekScalMat &edgemat = *edgeExp->GetLocMatrix(mkey);
......
......@@ -978,7 +978,7 @@ namespace Nektar
Array<OneD, int> sign;
StdRegions::VarCoeffMap varcoeffs;
varcoeffs[StdRegions::eVarCoeffPrimative] = primCoeffs;
varcoeffs[StdRegions::eVarCoeffMass] = primCoeffs;
StdRegions::ExpansionType expType =
faceExp->DetExpansionType();
......
......@@ -1297,7 +1297,8 @@ cout<<"deps/dx ="<<inarray_d0[i]<<" deps/dy="<<inarray_d1[i]<<endl;
{
case StdRegions::eMass:
{
if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
if((m_metricinfo->GetGtype() == SpatialDomains::eDeformed)||
(mkey.GetNVarCoeff()))
{
fac = 1.0;
goto UseLocRegionsMatrix;
......@@ -1311,7 +1312,8 @@ cout<<"deps/dx ="<<inarray_d0[i]<<" deps/dy="<<inarray_d1[i]<<endl;
break;
case StdRegions::eInvMass:
{
if(m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
if((m_metricinfo->GetGtype() == SpatialDomains::eDeformed)||
(mkey.GetNVarCoeff()))
{
NekDouble one = 1.0;
StdRegions::StdMatrixKey masskey(StdRegions::eMass,DetExpansionType(),
......
......@@ -1680,6 +1680,7 @@ namespace Nektar
int nblks = 2;
returnval = MemoryManager<DNekScalBlkMat>::AllocateSharedPtr(nblks, nblks, exp_size, exp_size); //Really need a constructor which takes Arrays
NekDouble factor = 1.0;
MatrixStorage AMatStorage = eFULL;
switch(mkey.GetMatrixType())
{
......@@ -1726,7 +1727,7 @@ namespace Nektar
NekDouble invfactor = 1.0/factor;
NekDouble one = 1.0;
DNekScalMat &mat = *GetLocMatrix(mkey);
DNekMatSharedPtr A = MemoryManager<DNekMat>::AllocateSharedPtr(nbdry,nbdry);
DNekMatSharedPtr A = MemoryManager<DNekMat>::AllocateSharedPtr(nbdry,nbdry,AMatStorage);
DNekMatSharedPtr B = MemoryManager<DNekMat>::AllocateSharedPtr(nbdry,nint);
DNekMatSharedPtr C = MemoryManager<DNekMat>::AllocateSharedPtr(nint,nbdry);
DNekMatSharedPtr D = MemoryManager<DNekMat>::AllocateSharedPtr(nint,nint);
......
......@@ -606,13 +606,65 @@ namespace Nektar
}
/****
* STEP 1.4: Check for singular system and add pinning Dirichlet vertex
*/
// Check between processes if the whole system is singular
int n = m_comm->GetSize();
int p = m_comm->GetRank();
int s = (systemSingular ? 1 : 0);
vCommRow->AllReduce(s, LibUtilities::ReduceMin);
systemSingular = (s == 1 ? true : false);
// Count the number of boundary regions on each process
Array<OneD, int> bccounts(n, 0);
bccounts[p] = bndCondExp.num_elements();
vCommRow->AllReduce(bccounts, LibUtilities::ReduceSum);
// Find the process rank with the maximum number of boundary regions
int maxBCIdx = Vmath::Imax(n, bccounts, 1);
// If the system is singular, the process with the maximum number of
// BCs will set a Dirichlet vertex to make system non-singular.
// Note: we find the process with maximum boundary regions to ensure
// we do not try to set a Dirichlet vertex on a partition with no
// intersection with the boundary.
if(systemSingular == true && checkIfSystemSingular && maxBCIdx == p)
{
if(m_session->DefinesParameter("SingularElement"))
{
int s_eid;
m_session->LoadParameter("SingularElement", s_eid);
ASSERTL1(s_eid < locExpVector.size(),"SingularElement Parameter is too large");
meshVertId = locExpVector[s_eid]->GetGeom2D()->GetVid(0);
}
else if (m_session->DefinesParameter("SingularVertex"))
{
m_session->LoadParameter("SingularVertex", meshVertId);
}
else
{
//last region i and j=0 edge
bndSegExp = boost::dynamic_pointer_cast<LocalRegions::SegExp>(bndCondExp[bndCondExp.num_elements()-1]->GetExp(0));
//first vertex 0 of the edge
meshVertId = (bndSegExp->GetGeom1D())->GetVid(0);
}
if(ReorderedGraphVertId[0].count(meshVertId) == 0)
{
ReorderedGraphVertId[0][meshVertId] = graphVertId++;
}
}
/**
* STEP 1.5: Exchange Dirichlet mesh vertices between processes and
* check for singular problems.
*/
// Collate information on Dirichlet vertices from all processes
int n = m_comm->GetSize();
int p = m_comm->GetRank();
Array<OneD, int> counts (n, 0);
Array<OneD, int> offsets(n, 0);
counts[p] = ReorderedGraphVertId[0].size();
......@@ -711,50 +763,6 @@ namespace Nektar
extraDirVerts.insert(vertids[i]);
}
// Check between processes if the whole system is singular
int s = (systemSingular ? 1 : 0);
vCommRow->AllReduce(s, LibUtilities::ReduceMin);
systemSingular = (s == 1 ? true : false);
// Count the number of boundary regions on each process
Array<OneD, int> bccounts(n, 0);
bccounts[p] = bndCondExp.num_elements();
vCommRow->AllReduce(bccounts, LibUtilities::ReduceSum);
// Find the process rank with the maximum number of boundary regions
int maxBCIdx = Vmath::Imax(n, bccounts, 1);
// If the system is singular, the process with the maximum number of
// BCs will set a Dirichlet vertex to make system non-singular.
// Note: we find the process with maximum boundary regions to ensure
// we do not try to set a Dirichlet vertex on a partition with no
// intersection with the boundary.
if(systemSingular == true && checkIfSystemSingular && maxBCIdx == p)
{
if(m_session->DefinesParameter("SingularElement"))
{
int s_eid;
m_session->LoadParameter("SingularElement", s_eid);
ASSERTL1(s_eid < locExpVector.size(),"SingularElement Parameter is too large");
meshVertId = locExpVector[s_eid]->GetGeom2D()->GetVid(0);
}
else
{
//last region i and j=0 edge
bndSegExp = boost::dynamic_pointer_cast<LocalRegions::SegExp>(bndCondExp[bndCondExp.num_elements()-1]->GetExp(0));
//first vertex 0 of the edge
meshVertId = (bndSegExp->GetGeom1D())->GetVid(0);
}
if(ReorderedGraphVertId[0].count(meshVertId) == 0)
{
ReorderedGraphVertId[0][meshVertId] = graphVertId++;
}
}
firstNonDirGraphVertId = graphVertId;
/**
......
......@@ -1097,7 +1097,7 @@ namespace Nektar
}
/**
* @brief Set up a list of elemeent IDs and edge IDs that link to the
* @brief Set up a list of element IDs and edge IDs that link to the
* boundary conditions.
*/
void DisContField2D::v_GetBoundaryToElmtMap(
......@@ -1690,7 +1690,7 @@ namespace Nektar
std::vector<SpatialDomains::
FieldDefinitionsSharedPtr> FieldDef;
std::vector<std::vector<NekDouble> > FieldData;
std::vector<std::vector<NekDouble> > FieldData;
m_graph->Import(filebcs,FieldDef, FieldData);
......
......@@ -1396,12 +1396,21 @@ namespace Nektar
SpatialDomains::RobinBoundaryCondition
>(m_bndConditions[i])->m_robinFunction;
LibUtilities::Equation coeff =
boost::static_pointer_cast<
SpatialDomains::RobinBoundaryCondition
>(m_bndConditions[i])->m_robinPrimitiveCoeff;
condition.Evaluate(x0,x1,x2,time,locExpList->UpdatePhys());
locExpList->IProductWRTBase(locExpList->GetPhys(),
locExpList->UpdateCoeffs());
/// \todo RobinPrimitiveCoeff forgotten? - PB
// put primitive coefficient into the physical space
// storage
coeff.Evaluate(x0,x1,x2,time,
locExpList->UpdatePhys());
}
else
{
......
......@@ -297,14 +297,14 @@ namespace Nektar
}
}
/**
/**
* Fills the list of local expansions with the segments in one
* subdomain specified in an inputfile by \a domain. This
* CompositeMap contains a list of Composites which define the
* subdomains.
* subdomain specified in an inputfile by \a domain. This
* CompositeMap contains a list of Composites which define the
* subdomains.
* @param domain A domain, comprising of one or more composite
* regions.
* @param i Index of currently processed subdomain
* @param i Index of currently processed subdomain
* @param graph1D A mesh, containing information about the
* domain and the spectral/hp element expansion.
* @param DeclareCoeffPhysArrays If true, create general segment expansions
......
......@@ -170,8 +170,7 @@ namespace Nektar
* Given a block matrix, construct a global matrix system according to
* a local to global mapping. #m_linSys is constructed by
* AssembleFullMatrix().
* @param mkey Associated linear system key.
* @param Mat Block matrix.
* @param pkey Associated linear system key.
* @param locToGloMap Local to global mapping.
*/
GlobalLinSys::GlobalLinSys(const GlobalLinSysKey &pKey,
......
......@@ -867,9 +867,11 @@ namespace Nektar
else
{
// Do matrix multiply locally
NekVector<NekDouble> loc(nLocal, m_wsp, eWrapper);
m_locToGloMap->GlobalToLocalBnd(pInput, m_wsp);
NekVector<NekDouble> loc(nLocal, m_wsp, eWrapper);
loc = (*m_schurCompl)*loc;
m_locToGloMap->AssembleBnd(m_wsp, pOutput);
}
}
......
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