Skip to content

GitLab

Commit 89dc3fef authored by Dave Moxey's avatar Dave Moxey
Browse files

Merge branch 'master' into feature/cfs-optimise

parents 8120ba23 966d08c6
Hide whitespace changes
Inline Side-by-side
library/LibUtilities/Interpreter/AnalyticExpressionEvaluator.cpp
View file @ 89dc3fef
......@@ -684,7 +684,9 @@ namespace Nektar
boost::algorithm::trim(valueStr);
const parser_id parserID = location->value.id();
#if defined(NEKTAR_DEBUG) || defined(NEKTAR_FULLDEBUG)
const int num_children = location->children.size();
#endif
if (parserID == AnalyticExpression::constantID)
{
......
library/LocalRegions/HexExp.cpp
View file @ 89dc3fef
......@@ -1289,7 +1289,28 @@ namespace Nektar
SpatialDomains::GeomType type = geomFactors->GetGtype();
const Array<TwoD, const NekDouble> & gmat = geomFactors->GetGmat();
const Array<OneD, const NekDouble> & jac = geomFactors->GetJac();
int nqe = m_base[0]->GetNumPoints()*m_base[1]->GetNumPoints();
int nqe0 = m_base[0]->GetNumPoints();
int nqe1 = m_base[1]->GetNumPoints();
int nqe2 = m_base[2]->GetNumPoints();
int nqe01 = nqe0*nqe1;
int nqe02 = nqe0*nqe2;
int nqe12 = nqe1*nqe2;
int nqe;
if (face == 0 || face == 5)
{
nqe = nqe01;
}
else if (face == 1 || face == 3)
{
nqe = nqe02;
}
else
{
nqe = nqe12;
}
int vCoordDim = GetCoordim();
m_faceNormals[face] = Array<OneD, Array<OneD, NekDouble> >(vCoordDim);
......@@ -1357,176 +1378,90 @@ namespace Nektar
Vmath::Smul(nqe,fac,normal[i],1,normal[i],1);
}
}
}
else // Set up deformed normals
{
int j, k;
int nquad0 = geomFactors->GetPointsKey(0).GetNumPoints();
int nquad1 = geomFactors->GetPointsKey(1).GetNumPoints();
int nquad2 = geomFactors->GetPointsKey(2).GetNumPoints();
int nqtot = nquad0*nquad1;
LibUtilities::PointsKey points0;
LibUtilities::PointsKey points1;
{
int j, k;
int nq;
Array<OneD,NekDouble> work(nqe,0.0);
Array<OneD,NekDouble> normals(vCoordDim*nquad0*nquad1,0.0);
//Array<OneD,NekDouble> normals(vCoordDim*nqe,0.0);
// Extract Jacobian along face and recover local
// derivates (dx/dr) for polynomial interpolation by
// multiplying m_gmat by jacobian
switch(face)
{
case 0:
for(j = 0; j < nquad0*nquad1; ++j)
{
normals[j] = -gmat[2][j]*jac[j];
normals[nqtot+j] = -gmat[5][j]*jac[j];
normals[2*nqtot+j] = -gmat[8][j]*jac[j];
}
points0 = geomFactors->GetPointsKey(0);
points1 = geomFactors->GetPointsKey(1);
nq=points0.GetNumPoints()*points1.GetNumPoints();
// interpolate Jacobian and invert
LibUtilities::Interp2D(points0,points1,jac,m_base[0]->GetPointsKey(),m_base[1]->GetPointsKey(),work);
Vmath::Sdiv(nq,1.0,&work[0],1,&work[0],1);
// interpolate
for(i = 0; i < GetCoordim(); ++i)
{
LibUtilities::Interp2D(points0,points1,&normals[i*nq],m_base[0]->GetPointsKey(),m_base[1]->GetPointsKey(),&normal[i][0]);
Vmath::Vmul(nqe,work,1,normal[i],1,normal[i],1);
}
break;
case 1:
for (j=0; j< nquad0; ++j)
{
for(k=0; k<nquad2; ++k)
{
case 0:
for(j = 0; j < nqe; ++j)
{
normals[j+k*nquad0] = -gmat[1][j+nquad0*nquad1*k]*jac[j+nquad0*nquad1*k];
normals[nqtot+j+k*nquad0] = -gmat[4][j+nquad0*nquad1*k]*jac[j+nquad0*nquad1*k];
normals[2*nqtot+j+k*nquad0] = -gmat[7][j+nquad0*nquad1*k]*jac[j+nquad0*nquad1*k];
}
}
points0 = geomFactors->GetPointsKey(0);
points1 = geomFactors->GetPointsKey(2);
nq=points0.GetNumPoints()*points1.GetNumPoints();
// interpolate Jacobian and invert
LibUtilities::Interp2D(points0,points1,jac,m_base[0]->GetPointsKey(),m_base[2]->GetPointsKey(),work);
Vmath::Sdiv(nq,1.0,&work[0],1,&work[0],1);
// interpolate
for(i = 0; i < GetCoordim(); ++i)
{
LibUtilities::Interp2D(points0,points1,&normals[i*nq],m_base[0]->GetPointsKey(),m_base[2]->GetPointsKey(),&normal[i][0]);
Vmath::Vmul(nqe,work,1,normal[i],1,normal[i],1);
}
break;
case 2:
for (j=0; j< nquad1; ++j)
{
for(k=0; k<nquad2; ++k)
normal[0][j] = -gmat[2][j]*jac[j];
normal[1][j] = -gmat[5][j]*jac[j];
normal[2][j] = -gmat[8][j]*jac[j];
}
break;
case 1:
for (j = 0; j< nqe0; ++j)
{
normals[j+k*nquad0] = gmat[0][nquad0-1+nquad0*j+nquad0*nquad1*k]*jac[nquad0-1+nquad0*j+nquad0*nquad1*k];
normals[nqtot+j+k*nquad0] = gmat[3][nquad0-1+nquad0*j+nquad0*nquad1*k]*jac[nquad0-1+nquad0*j+nquad0*nquad1*k];
normals[2*nqtot+j+k*nquad0] = gmat[6][nquad0-1+nquad0*j+nquad0*nquad1*k]*jac[nquad0-1+nquad0*j+nquad0*nquad1*k];
}
}
points0 = geomFactors->GetPointsKey(1);
points1 = geomFactors->GetPointsKey(2);
nq=points0.GetNumPoints()*points1.GetNumPoints();
// interpolate Jacobian and invert
LibUtilities::Interp2D(points0,points1,jac,m_base[1]->GetPointsKey(),m_base[2]->GetPointsKey(),work);
Vmath::Sdiv(nq,1.0,&work[0],1,&work[0],1);
// interpolate
for(i = 0; i < GetCoordim(); ++i)
{
LibUtilities::Interp2D(points0,points1,&normals[i*nq],m_base[1]->GetPointsKey(),m_base[2]->GetPointsKey(),&normal[i][0]);
Vmath::Vmul(nqe,work,1,normal[i],1,normal[i],1);
}
break;
case 3:
for (j=0; j< nquad0; ++j)
{
for(k=0; k<nquad2; ++k)
for(k = 0; k < nqe2; ++k)
{
normal[0][j+k*nqe0] = -gmat[1][j+nqe01*k]*jac[j+nqe01*k];
normal[1][j+k*nqe0] = -gmat[4][j+nqe01*k]*jac[j+nqe01*k];
normal[2][j+k*nqe0] = -gmat[7][j+nqe01*k]*jac[j+nqe01*k];
}
}
break;
case 2:
for (j=0; j< nqe1; ++j)
{
normals[j+k*nquad0] = gmat[1][nquad0*(nquad1-1)+j+nquad0*nquad1*k]*jac[nquad0*(nquad1-1)+j+nquad0*nquad1*k];
normals[nqtot+j+k*nquad0] = gmat[4][nquad0*(nquad1-1)+j+nquad0*nquad1*k]*jac[nquad0*(nquad1-1)+j+nquad0*nquad1*k];
normals[2*nqtot+j+k*nquad0] = gmat[7][nquad0*(nquad1-1)+j+nquad0*nquad1*k]*jac[nquad0*(nquad1-1)+j+nquad0*nquad1*k];
}
}
points0 = geomFactors->GetPointsKey(0);
points1 = geomFactors->GetPointsKey(2);
nq=points0.GetNumPoints()*points1.GetNumPoints();
// interpolate Jacobian and invert
LibUtilities::Interp2D(points0,points1,jac,m_base[0]->GetPointsKey(),m_base[2]->GetPointsKey(),work);
Vmath::Sdiv(nq,1.0,&work[0],1,&work[0],1);
// interpolate
for(i = 0; i < GetCoordim(); ++i)
{
LibUtilities::Interp2D(points0,points1,&normals[i*nq],m_base[0]->GetPointsKey(),m_base[2]->GetPointsKey(),&normal[i][0]);
Vmath::Vmul(nqe,work,1,normal[i],1,normal[i],1);
}
break;
case 4:
for (j=0; j< nquad0; ++j)
{
for(k=0; k<nquad2; ++k)
for(k=0; k<nqe2; ++k)
{
normal[0][j+k*nqe0] = gmat[0][nqe0-1+nqe0*j+nqe01*k]*jac[nqe0-1+nqe0*j+nqe01*k];
normal[1][j+k*nqe0] = gmat[3][nqe0-1+nqe0*j+nqe01*k]*jac[nqe0-1+nqe0*j+nqe01*k];
normal[2][j+k*nqe0] = gmat[6][nqe0-1+nqe0*j+nqe01*k]*jac[nqe0-1+nqe0*j+nqe01*k];
}
}
break;
case 3:
for (j=0; j< nqe0; ++j)
{
normals[j+k*nquad0] = -gmat[0][j*nquad0+nquad0*nquad1*k]*jac[j*nquad0+nquad0*nquad1*k];
normals[nqtot+j+k*nquad0] = -gmat[3][j*nquad0+nquad0*nquad1*k]*jac[j*nquad0+nquad0*nquad1*k];
normals[2*nqtot+j+k*nquad0] = -gmat[6][j*nquad0+nquad0*nquad1*k]*jac[j*nquad0+nquad0*nquad1*k];
}
}
points0 = geomFactors->GetPointsKey(1);
points1 = geomFactors->GetPointsKey(2);
nq=points0.GetNumPoints()*points1.GetNumPoints();
// interpolate Jacobian and invert
LibUtilities::Interp2D(points0,points1,jac,m_base[1]->GetPointsKey(),m_base[2]->GetPointsKey(),work);
Vmath::Sdiv(nq,1.0,&work[0],1,&work[0],1);
for(k=0; k<nqe2; ++k)
{
normal[0][j+k*nqe0] = gmat[1][nqe0*(nqe1-1)+j+nqe0*nqe1*k]*jac[nqe0*(nqe1-1)+j+nqe01*k];
normal[1][j+k*nqe0] = gmat[4][nqe0*(nqe1-1)+j+nqe0*nqe1*k]*jac[nqe0*(nqe1-1)+j+nqe01*k];
normal[1][j+k*nqe0] = gmat[7][nqe0*(nqe1-1)+j+nqe0*nqe1*k]*jac[nqe0*(nqe1-1)+j+nqe01*k];
}
}
break;
case 4:
for (j=0; j< nqe0; ++j)
{
for(k=0; k<nqe2; ++k)
{
normal[0][j+k*nqe0] = -gmat[0][j*nqe0+nqe01*k]*jac[j*nqe0+nqe01*k];
normal[1][j+k*nqe0] = -gmat[3][j*nqe0+nqe01*k]*jac[j*nqe0+nqe01*k];
normal[2][j+k*nqe0] = -gmat[6][j*nqe0+nqe01*k]*jac[j*nqe0+nqe01*k];
}
}
break;
case 5:
for (j=0; j< nqe01; ++j)
{
normal[0][j] = gmat[2][j+nqe01*(nqe2-1)]*jac[j+nqe01*(nqe2-1)];
normal[1][j] = gmat[5][j+nqe01*(nqe2-1)]*jac[j+nqe01*(nqe2-1)];
normal[2][j] = gmat[8][j+nqe01*(nqe2-1)]*jac[j+nqe01*(nqe2-1)];
}
break;
default:
ASSERTL0(false,"face is out of range (face < 5)");
}
// interpolate
for(i = 0; i < GetCoordim(); ++i)
{
LibUtilities::Interp2D(points0,points1,&normals[i*nq],m_base[1]->GetPointsKey(),m_base[2]->GetPointsKey(),&normal[i][0]);
Vmath::Vmul(nqe,work,1,normal[i],1,normal[i],1);
}
break;
case 5:
for (j=0; j< nquad0*nquad1; ++j)
{
normals[j] = gmat[2][j+nquad0*nquad1*(nquad2-1)]*jac[j+nquad0*nquad1*(nquad2-1)];
normals[nqtot+j] = gmat[5][j+nquad0*nquad1*(nquad2-1)]*jac[j+nquad0*nquad1*(nquad2-1)];
normals[2*nqtot+j] = gmat[8][j+nquad0*nquad1*(nquad2-1)]*jac[j+nquad0*nquad1*(nquad2-1)];
}
points0 = geomFactors->GetPointsKey(0);
points1 = geomFactors->GetPointsKey(1);
nq=points0.GetNumPoints()*points1.GetNumPoints();
// interpolate Jacobian and invert
LibUtilities::Interp2D(points0,points1,jac,m_base[0]->GetPointsKey(),m_base[1]->GetPointsKey(),work);
Vmath::Sdiv(nq,1.0,&work[0],1,&work[0],1);
Vmath::Sdiv(nqe,1.0,&jac[0],1,&work[0],1);
// interpolate
for(i = 0; i < GetCoordim(); ++i)
{
LibUtilities::Interp2D(points0,points1,&normals[i*nq],m_base[0]->GetPointsKey(),m_base[1]->GetPointsKey(),&normal[i][0]);
Vmath::Vmul(nqe,work,1,normal[i],1,normal[i],1);
}
break;
default:
ASSERTL0(false,"face is out of range (face < 5)");
}
// interpolate
for(i = 0; i < GetCoordim(); ++i)
{
Vmath::Vmul(nqe,work,1,normal[i],1,normal[i],1);
}
//normalise normal vectors
Vmath::Zero(nqe,work,1);
......@@ -1541,7 +1476,7 @@ namespace Nektar
for(i = 0; i < GetCoordim(); ++i)
{
Vmath::Vmul(nqe,normal[i],1,work,1,normal[i],1);
}
}
}
}
......
library/MultiRegions/GlobalLinSysIterative.cpp
View file @ 89dc3fef
......@@ -141,6 +141,10 @@ namespace Nektar
const AssemblyMapSharedPtr &plocToGloMap,
const int nDir)
{
// Get the communicator for performing data exchanges
LibUtilities::CommSharedPtr vComm
= m_expList.lock()->GetComm()->GetRowComm();
// Get vector sizes
int nNonDir = nGlobal - nDir;
Array<OneD, NekDouble> tmp;
......@@ -161,9 +165,12 @@ namespace Nektar
// check the input vector (rhs) is not zero
NekDouble rhsNorm = Vmath::Dot(nNonDir,
pInput + nDir,
pInput + nDir);
NekDouble rhsNorm = Vmath::Dot2(nNonDir,
pInput + nDir,
pInput + nDir,
m_map + nDir);
vComm->AllReduce(rhsNorm, Nektar::LibUtilities::ReduceSum);
if (rhsNorm < NekConstants::kNekZeroTol)
{
......@@ -188,19 +195,25 @@ namespace Nektar
// \alpha_i = \tilda{x_i}^T b^n
// projected x, px = \sum \alpha_i \tilda{x_i}
Array<OneD, NekDouble> alpha (m_prevLinSol.size(), 0.0);
for (int i = 0; i < m_prevLinSol.size(); i++)
{
NekDouble alphai = Vmath::Dot(nNonDir,
m_prevLinSol[i],
pInput + nDir);
alpha[i] = Vmath::Dot2(nNonDir,
m_prevLinSol[i],
pInput + nDir,
m_map + nDir);
}
vComm->AllReduce(alpha, Nektar::LibUtilities::ReduceSum);
if (alphai < NekConstants::kNekZeroTol)
for (int i = 0; i < m_prevLinSol.size(); i++)
{
if (alpha[i] < NekConstants::kNekZeroTol)
{
continue;
}
NekVector<NekDouble> xi (nNonDir, m_prevLinSol[i], eWrapper);
px += alphai * xi;
px += alpha[i] * xi;
}
// pb = b^n - A px
......@@ -235,6 +248,10 @@ namespace Nektar
const Array<OneD,const NekDouble> &in,
const int nDir)
{
// Get the communicator for performing data exchanges
LibUtilities::CommSharedPtr vComm
= m_expList.lock()->GetComm()->GetRowComm();
// Get vector sizes
int nNonDir = nGlobal - nDir;
......@@ -243,9 +260,11 @@ namespace Nektar
v_DoMatrixMultiply(in, tmpAx_s);
const NekDouble anorm_sq = Vmath::Dot(nNonDir,
in + nDir,
tmpAx_s + nDir);
NekDouble anorm_sq = Vmath::Dot2(nNonDir,
in + nDir,
tmpAx_s + nDir,
m_map + nDir);
vComm->AllReduce(anorm_sq, Nektar::LibUtilities::ReduceSum);
return std::sqrt(anorm_sq);
}
......@@ -261,11 +280,16 @@ namespace Nektar
// Get vector sizes
int nNonDir = nGlobal - nDir;
// Get the communicator for performing data exchanges
LibUtilities::CommSharedPtr vComm
= m_expList.lock()->GetComm()->GetRowComm();
// Check the solution is non-zero
NekDouble solNorm = Vmath::Dot(nNonDir,
newX + nDir,
newX + nDir);
NekDouble solNorm = Vmath::Dot2(nNonDir,
newX + nDir,
newX + nDir,
m_map + nDir);
vComm->AllReduce(solNorm, Nektar::LibUtilities::ReduceSum);
if (solNorm < NekConstants::kNekZeroTol)
{
......@@ -293,19 +317,26 @@ namespace Nektar
{
v_DoMatrixMultiply(newX, tmpAx_s);
}
Array<OneD, NekDouble> alpha (m_prevLinSol.size(), 0.0);
for (int i = 0; i < m_prevLinSol.size(); i++)
{
NekDouble alphai = Vmath::Dot(nNonDir,
m_prevLinSol[i],
tmpAx_s + nDir);
alpha[i] = Vmath::Dot2(nNonDir,
m_prevLinSol[i],
tmpAx_s + nDir,
m_map + nDir);
}
vComm->AllReduce(alpha, Nektar::LibUtilities::ReduceSum);
if (alphai < NekConstants::kNekZeroTol)
for (int i = 0; i < m_prevLinSol.size(); i++)
{
if (alpha[i] < NekConstants::kNekZeroTol)
{
continue;
}
NekVector<NekDouble> xi (nNonDir, m_prevLinSol[i], eWrapper);
px -= alphai * xi;
px -= alpha[i] * xi;
}
......
solvers/CardiacEPSolver/EquationSystems/Monodomain.cpp
View file @ 89dc3fef
......@@ -109,24 +109,39 @@ namespace Nektar
int nq = m_fields[0]->GetNpoints();
Array<OneD, NekDouble> vTemp;
// Allocate storage for variable coeffs and initialize to 1.
for (int i = 0; i < m_spacedim; ++i)
{
m_vardiff[varCoeffEnum[i]] = Array<OneD, NekDouble>(nq, 1.0);
}
// Apply intensity map (range d_min -> d_max)
if (m_session->DefinesFunction("IsotropicConductivity"))
{
if (m_session->DefinesCmdLineArgument("verbose"))
{
cout << "Loading Isotropic Conductivity map." << endl;
}
EvaluateFunction(varName, vTemp, "IsotropicConductivity");
for (int i = 0; i < m_spacedim; ++i)
{
m_vardiff[varCoeffEnum[i]] = Array<OneD, NekDouble>(nq);
Vmath::Vcopy(nq, vTemp, 1, m_vardiff[varCoeffEnum[i]], 1);
Vmath::Vmul(nq, vTemp, 1, m_vardiff[varCoeffEnum[i]], 1, m_vardiff[varCoeffEnum[i]], 1);
}
}
else if (m_session->DefinesFunction(varCoeffs[0]))
// Apply fibre map (range 0 -> 1)
if (m_session->DefinesFunction(varCoeffs[0]))
{
if (m_session->DefinesCmdLineArgument("verbose"))
{
cout << "Loading Anisotropic Fibre map." << endl;
}
for (int i = 0; i < m_spacedim; ++i)
{
ASSERTL0(m_session->DefinesFunction(varCoeffs[i], varName),
"Function '" + varCoeffs[i] + "' not correctly defined.");
EvaluateFunction(varName, vTemp, varCoeffs[i]);
m_vardiff[varCoeffEnum[i]] = Array<OneD, NekDouble>(nq);
Vmath::Vcopy(nq, vTemp, 1, m_vardiff[varCoeffEnum[i]], 1);
Vmath::Vmul(nq, vTemp, 1, m_vardiff[varCoeffEnum[i]], 1, m_vardiff[varCoeffEnum[i]], 1);
}
}
......
solvers/IncNavierStokesSolver/CMakeLists.txt
View file @ 89dc3fef
......@@ -78,6 +78,7 @@ IF( NEKTAR_SOLVER_INCNAVIERSTOKES )
ADD_NEKTAR_TEST(Channel_Flow_3modes_Parallel)
ADD_NEKTAR_TEST(ChanFlow_m3_par)
ADD_NEKTAR_TEST(Hex_channel_m8_par)
ADD_NEKTAR_TEST(Hex_channel_m8_srhs_par)
ADD_NEKTAR_TEST(Tet_channel_m8_par)
ADD_NEKTAR_TEST(Tet_channel_m8_iter_ml_par)
ENDIF (NEKTAR_USE_MPI)
......
solvers/IncNavierStokesSolver/Tests/Hex_channel_m8_srhs_par.tst 0 → 100644
View file @ 89dc3fef
<?xml version="1.0" encoding="utf-8"?>
<test>
<description>3D channel flow, Hexahedral elements, P=8, Successive RHS(5), par(2)</description>
<executable>IncNavierStokesSolver</executable>
<parameters>Hex_channel_m8_srhs.xml</parameters>
<processes>2</processes>
<files>
<file description="Session File">Hex_channel_m8_srhs.xml</file>
</files>
<metrics>
<metric type="L2" id="1">
<value variable="u" tolerance="1e-8">2.73101e-11</value>
<value variable="v" tolerance="1e-8">2.23013e-11</value>
<value variable="w" tolerance="1e-8">6.29721e-10</value>
<value variable="p" tolerance="1e-8">1.94395e-09</value>
</metric>
<metric type="Linf" id="2">
<value variable="u" tolerance="1e-8">1.55008e-10</value>
<value variable="v" tolerance="1e-8">1.58963e-10</value>
<value variable="w" tolerance="1e-8">8.27366e-08</value>
<value variable="p" tolerance="1e-8">2.42343e-08</value>
</metric>
</metrics>
</test>
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