Commit ab244358 authored by Michael Turner's avatar Michael Turner
Browse files

Merge branch 'feature/jan_masterupdates' into 'master'

feature/varopti_updates

New node-colouring and streamlined GetFunctional()

See merge request !738
parents 767ff995 9f598de9
......@@ -93,6 +93,12 @@ v4.4.0
- Change variable names in mcf file to make more sense (!736)
- Fix issues in varopti module so that in can be compiled without meshgen on
(!736)
- Replace LAPACK Eigenvalue calculation with handwritten function in
varopti (!738)
- Improved node-colouring algorithm for better load-balancing
in varopti (!738)
- Simplified calculation of the energy functional in varopti for improved
performance (!738)
**FieldConvert:**
- Move all modules to a new library, FieldUtils, to support post-processing
......
......@@ -148,7 +148,8 @@ IF(NEKTAR_USE_MESHGEN)
ADD_NEKTAR_TEST (MeshGen/cylinder)
ADD_NEKTAR_TEST (MeshGen/sphere)
ADD_NEKTAR_TEST (MeshGen/2d-cad)
ADD_NEKTAR_TEST (MeshGen/2d-naca)
# varopti tests
ADD_NEKTAR_TEST_LENGTHY (MeshGen/varopti_cubesphere)
ADD_NEKTAR_TEST (MeshGen/2d-naca)
ENDIF()
# varopti tests
ADD_NEKTAR_TEST_LENGTHY (MeshGen/varopti_cubesphere)
......@@ -36,6 +36,7 @@
#ifndef UTILITIES_NEKMESH_NODEOPTI_HESSIAN
#define UTILITIES_NEKMESH_NODEOPTI_HESSIAN
namespace Nektar
{
namespace Utilities
......@@ -151,87 +152,107 @@ template <int DIM> void NodeOpti::MinEigen(NekDouble &val)
template <> void NodeOpti::MinEigen<2>(NekDouble &val)
{
Array<OneD, NekDouble> eigR(2);
Array<OneD, NekDouble> eigI(2);
NekMatrix<NekDouble> H(2, 2);
H(0, 0) = m_grad[2];
H(1, 0) = m_grad[3];
H(0, 1) = H(1, 0);
H(1, 1) = m_grad[4];
NekDouble H[2][2];
H[0][0] = m_grad[2];
H[1][0] = m_grad[3];
//H[0][1] = H[1][0];
H[1][1] = m_grad[4];
int nVel = 2;
char jobvl = 'N', jobvr = 'V';
int worklen = 8 * nVel, info;
//double eval[2]; // the eigenvalues
DNekMat eval(nVel, nVel, 0.0, eDIAGONAL);
DNekMat evec(nVel, nVel, 0.0, eFULL);
Array<OneD, NekDouble> vl(nVel * nVel);
Array<OneD, NekDouble> work(worklen);
Array<OneD, NekDouble> wi(nVel);
NekDouble D = (H[0][0] - H[1][1]) * (H[0][0] - H[1][1]) + 4.0 * H[1][0] * H[1][0];
NekDouble Dsqrt = sqrt(D);
Lapack::Dgeev(jobvl, jobvr, nVel, H.GetRawPtr(), nVel, &(eval.GetPtr())[0],
&wi[0], &vl[0], nVel, &(evec.GetPtr())[0], nVel, &work[0],
worklen, info);
ASSERTL0(!info, "dgeev failed");
int minI;
NekDouble tmp = std::numeric_limits<double>::max();
for (int i = 0; i < 2; i++)
{
if (eval(i, i) < tmp)
{
minI = i;
tmp = eval(i, i);
}
}
val = eval(minI, minI);
//eval[0] = (H[0][0] + H[1][1] + Dsqrt ) / 2.0;
val = (H[0][0] + H[1][1] - Dsqrt ) / 2.0; // the minimum Eigenvalue
}
template <> void NodeOpti::MinEigen<3>(NekDouble &val)
{
Array<OneD, NekDouble> eigR(3);
Array<OneD, NekDouble> eigI(3);
NekMatrix<NekDouble> H(3, 3);
H(0, 0) = m_grad[3];
H(1, 0) = m_grad[4];
H(0, 1) = H(1, 0);
H(2, 0) = m_grad[5];
H(0, 2) = H(2, 0);
H(1, 1) = m_grad[6];
H(2, 1) = m_grad[7];
H(1, 2) = H(2, 1);
H(2, 2) = m_grad[8];
int nVel = 3;
char jobvl = 'N', jobvr = 'V';
int worklen = 8 * nVel, info;
DNekMat eval(nVel, nVel, 0.0, eDIAGONAL);
DNekMat evec(nVel, nVel, 0.0, eFULL);
Array<OneD, NekDouble> vl(nVel * nVel);
Array<OneD, NekDouble> work(worklen);
Array<OneD, NekDouble> wi(nVel);
Lapack::Dgeev(jobvl, jobvr, nVel, H.GetRawPtr(), nVel, &(eval.GetPtr())[0],
&wi[0], &vl[0], nVel, &(evec.GetPtr())[0], nVel, &work[0],
worklen, info);
ASSERTL0(!info, "dgeev failed");
int minI;
NekDouble tmp = std::numeric_limits<double>::max();
for (int i = 0; i < 3; i++)
{
if (eval(i, i) < tmp)
NekDouble H[3][3];
H[0][0] = m_grad[3];
H[1][0] = m_grad[4];
H[0][1] = H[1][0];
H[2][0] = m_grad[5];
H[0][2] = H[2][0];
H[1][1] = m_grad[6];
H[2][1] = m_grad[7];
H[1][2] = H[2][1];
H[2][2] = m_grad[8];
//double eval[3]; // the eigenvalues
NekDouble p1 = H[0][1] * H[0][1] + H[0][2] * H[0][2] + H[1][2] * H[1][2];
if (p1 == 0.0) // H is diagonal
{
// find the minimum Eigenvalue
if(H[0][0] < H[1][1])
{
minI = i;
tmp = eval(i, i);
if(H[0][0] < H[2][2])
{
val = H[0][0];
}
else
{
val = H[2][2];
}
}
else
{
if(H[1][1] < H[2][2])
{
val = H[1][1];
}
else
{
val = H[2][2];
}
}
}
else
{
NekDouble q = (H[0][0] + H[1][1] + H[2][2]) / 3.0;
NekDouble p2 = (H[0][0] - q)*(H[0][0] - q)
+ (H[1][1] - q)*(H[1][1] - q)
+ (H[2][2] - q)*(H[2][2] - q)
+ 2.0 * p1;
NekDouble p = sqrt(p2 / 6.0);
NekDouble B[3][3]; // B = (1.0 / p) * (H - q * I) with I being the identity matrix
NekDouble pinv = 1.0 / p;
B[0][0] = pinv * (H[0][0] - q);
B[1][1] = pinv * (H[1][1] - q);
B[2][2] = pinv * (H[2][2] - q);
B[0][1] = pinv * H[0][1];
B[1][0] = B[0][1];
B[0][2] = pinv * H[0][2];
B[2][0] = B[0][2];
B[1][2] = pinv * H[1][2];
B[2][1] = B[1][2];
NekDouble r = Determinant<3>(B) / 2.0;
// In exact arithmetic for h symmetric matrix -1 <= r <= 1
// but computation error can leave it slightly outside this range.
NekDouble phi;
if (r <= -1)
{
phi = M_PI / 3.0;
}
else if (r >= 1)
{
phi = 0.0;
}
else
{
phi = acos(r) / 3.0;
}
val = eval(minI, minI);
// the eigenvalues satisfy eval[2] <= eval[1] <= eval[0]
//eval[0] = q + 2.0 * p * cos(phi);
val = q + 2.0 * p * cos(phi + (2.0*M_PI/3.0));
//eval[1] = 3.0 * q - eval[0] - eval[2]; // since trace(H) = eval[0] + eval[1] + eval[2]
}
}
}
......
......@@ -166,9 +166,28 @@ map<LibUtilities::ShapeType, DerivUtilSharedPtr> ProcessVarOpti::BuildDerivUtil(
return ret;
}
struct NodeComparator
{
const vector<int> & value_vector;
NodeComparator(const vector<int> & val_vec):
value_vector(val_vec) {}
bool operator()(int i1, int i2)
{
return value_vector[i1] > value_vector[i2];
}
};
vector<vector<NodeSharedPtr> > ProcessVarOpti::GetColouredNodes(
vector<ElUtilSharedPtr> elLock)
{
// create set of nodes to be ignored and hence not included in the coloursets
NodeSet ignoredNodes;
for (int i = 0; i < elLock.size(); i++)
{
......@@ -180,8 +199,7 @@ vector<vector<NodeSharedPtr> > ProcessVarOpti::GetColouredNodes(
}
}
// this figures out the dirclet nodes and colors the others into paralell
// sets
// create set of nodes which are at the boundary and hence not included in the colourset
NodeSet boundaryNodes;
switch (m_mesh->m_spaceDim)
......@@ -255,11 +273,16 @@ vector<vector<NodeSharedPtr> > ProcessVarOpti::GetColouredNodes(
}
default:
ASSERTL0(false,"space dim issue");
}
}
vector<NodeSharedPtr> remain;
//create vector of free nodes which "remain", hence will be included in the coloursets
vector<NodeSharedPtr> remainEdgeVertex;
vector<NodeSharedPtr> remainFace;
vector<NodeSharedPtr> remainVolume;
m_res->nDoF = 0;
// check if vertex nodes are in boundary or ignored nodes, otherwise add to EDGE-VERTEX remain nodes
NodeSet::iterator nit;
for (nit = m_mesh->m_vertexSet.begin(); nit != m_mesh->m_vertexSet.end();
++nit)
......@@ -268,7 +291,7 @@ vector<vector<NodeSharedPtr> > ProcessVarOpti::GetColouredNodes(
NodeSet::iterator nit3 = ignoredNodes.find(*nit);
if (nit2 == boundaryNodes.end() && nit3 == ignoredNodes.end())
{
remain.push_back(*nit);
remainEdgeVertex.push_back(*nit);
if ((*nit)->GetNumCadCurve() == 1)
{
m_res->nDoF++;
......@@ -284,6 +307,7 @@ vector<vector<NodeSharedPtr> > ProcessVarOpti::GetColouredNodes(
}
}
// check if edge nodes are in boundary or ignored nodes, otherwise add to EDGE-VERTEX remain nodes
EdgeSet::iterator eit;
for (eit = m_mesh->m_edgeSet.begin(); eit != m_mesh->m_edgeSet.end(); eit++)
{
......@@ -294,7 +318,7 @@ vector<vector<NodeSharedPtr> > ProcessVarOpti::GetColouredNodes(
NodeSet::iterator nit3 = ignoredNodes.find(n[j]);
if (nit2 == boundaryNodes.end() && nit3 == ignoredNodes.end())
{
remain.push_back(n[j]);
remainEdgeVertex.push_back(n[j]);
if (n[j]->GetNumCadCurve() == 1)
{
m_res->nDoF++;
......@@ -311,6 +335,7 @@ vector<vector<NodeSharedPtr> > ProcessVarOpti::GetColouredNodes(
}
}
// check if face nodes are in boundary or ignored nodes, otherwise add to FACE remain nodes
FaceSet::iterator fit;
for (fit = m_mesh->m_faceSet.begin(); fit != m_mesh->m_faceSet.end(); fit++)
{
......@@ -320,7 +345,7 @@ vector<vector<NodeSharedPtr> > ProcessVarOpti::GetColouredNodes(
NodeSet::iterator nit3 = ignoredNodes.find((*fit)->m_faceNodes[j]);
if (nit2 == boundaryNodes.end() && nit3 == ignoredNodes.end())
{
remain.push_back((*fit)->m_faceNodes[j]);
remainFace.push_back((*fit)->m_faceNodes[j]);
if ((*fit)->m_faceNodes[j]->GetNumCADSurf() == 1)
{
m_res->nDoF += 2;
......@@ -333,6 +358,7 @@ vector<vector<NodeSharedPtr> > ProcessVarOpti::GetColouredNodes(
}
}
// check if volume nodes are in boundary or ignored nodes, otherwise add to VOLUME remain nodes
for (int i = 0; i < m_mesh->m_element[m_mesh->m_expDim].size(); i++)
{
vector<NodeSharedPtr> ns =
......@@ -343,49 +369,129 @@ vector<vector<NodeSharedPtr> > ProcessVarOpti::GetColouredNodes(
NodeSet::iterator nit3 = ignoredNodes.find(ns[j]);
if (nit2 == boundaryNodes.end() && nit3 == ignoredNodes.end())
{
remain.push_back(ns[j]);
remainVolume.push_back(ns[j]);
m_res->nDoF += m_mesh->m_spaceDim;
}
}
}
m_res->n = remain.size();
// size of all free nodes to be included in the coloursets
m_res->n = remainEdgeVertex.size()
+ remainFace.size() + remainVolume.size();
// data structure for coloursets, that will ultimately contain all free nodes
vector<vector<NodeSharedPtr> > ret;
vector<vector<NodeSharedPtr> > retPart;
// edge and vertex nodes
// create vector num_el of number of associated elements of each node
vector<int> num_el(remainEdgeVertex.size());
for (int i = 0; i < remainEdgeVertex.size(); i++)
{
//try to find node within all elements
NodeElMap::iterator it = m_nodeElMap.find(remainEdgeVertex[i]->m_id);
vector<ElUtilSharedPtr> &elUtils = it->second;
num_el[i] = elUtils.size();
}
// finding the permutation according to num_el
vector<int> permNode(remainEdgeVertex.size());
for (int i = 0; i < remainEdgeVertex.size(); ++i)
{
permNode[i] = i;
}
std::sort(permNode.begin(), permNode.end(), NodeComparator(num_el));
// applying the permutation to remainEdgeVertex
vector<NodeSharedPtr> remainEdgeVertexSort(remainEdgeVertex.size());
for (int i = 0; i < remainEdgeVertex.size(); ++i)
{
int j = permNode[i];
remainEdgeVertexSort[i] = remainEdgeVertex[j];
}
retPart = CreateColoursets(remainEdgeVertexSort);
if(m_mesh->m_verbose)
{
printf("\nNumber of Edge/Vertex Coloursets: %i\n", retPart.size());
}
for (int i = 0; i < retPart.size(); i++)
{
ret.push_back(retPart[i]);
}
// face nodes
retPart = CreateColoursets(remainFace);
if(m_mesh->m_verbose)
{
printf("\nNumber of Face Coloursets: %i\n" ,retPart.size());
}
for (int i = 0; i < retPart.size(); i++)
{
ret.push_back(retPart[i]);
}
// volume nodes
retPart = CreateColoursets(remainVolume);
if(m_mesh->m_verbose)
{
printf("\nNumber of Volume Coloursets: %i\n", retPart.size());
}
for (int i = 0; i < retPart.size(); i++)
{
ret.push_back(retPart[i]);
}
if(m_mesh->m_verbose)
{
cout << endl;
}
return ret;
}
vector<vector<NodeSharedPtr> > ProcessVarOpti::CreateColoursets(
vector<NodeSharedPtr> remain)
{
vector<vector<NodeSharedPtr> > retPart;
// loop until all free nodes have been sorted
while (remain.size() > 0)
{
vector<NodeSharedPtr> layer;
vector<NodeSharedPtr> layer; // one colourset
set<int> locked;
set<int> completed;
for (int i = 0; i < remain.size(); i++)
{
NodeElMap::iterator it = m_nodeElMap.find(remain[i]->m_id);
ASSERTL0(it != m_nodeElMap.end(), "could not find");
NodeElMap::iterator it = m_nodeElMap.find(remain[i]->m_id); //try to find node within all elements
ASSERTL0(it != m_nodeElMap.end(), "could not find node");
vector<ElUtilSharedPtr> &elUtils = it->second;
vector<ElUtilSharedPtr> &elUtils = it->second; // identify the vector of all associated elements of the node
bool islocked = false;
for (int j = 0; j < elUtils.size(); j++)
bool islocked = false; // suppose node is not locked
for (int j = 0; j < elUtils.size(); j++) // loop over all associated elements of the node
{
set<int>::iterator sit = locked.find(elUtils[j]->GetId());
if (sit != locked.end())
set<int>::iterator sit = locked.find(elUtils[j]->GetId()); // check all nodes of the element
if (sit != locked.end()) //if node is within the set of locked nodes
{
islocked = true;
break;
islocked = true; // if yes, flag node as locked
break; // and go to next node
}
}
if (!islocked)
if (!islocked) // if the node is not locked
{
layer.push_back(remain[i]);
completed.insert(remain[i]->m_id);
for (int j = 0; j < elUtils.size(); j++)
layer.push_back(remain[i]); // insert node into colourset
completed.insert(remain[i]->m_id); // insert sorted node into "completed" list
for (int j = 0; j < elUtils.size(); j++) // loop over all other nodes of the same element
{
locked.insert(elUtils[j]->GetId());
locked.insert(elUtils[j]->GetId()); // and flag these nodes as locked
}
}
}
// identify nodes which are not sorted, yet and create new "remain" vector
vector<NodeSharedPtr> tmp = remain;
remain.clear();
for (int i = 0; i < tmp.size(); i++)
......@@ -396,22 +502,21 @@ vector<vector<NodeSharedPtr> > ProcessVarOpti::GetColouredNodes(
remain.push_back(tmp[i]);
}
}
ret.push_back(layer);
// include layer or colourset into vector of coloursets
retPart.push_back(layer);
// print out progress
if(m_mesh->m_verbose)
{
LibUtilities::PrintProgressbar(m_res->n - remain.size(), m_res->n, "Node Coloring");
}
}
if(m_mesh->m_verbose)
{
cout << endl;
}
return ret;
return retPart;
}
void ProcessVarOpti::GetElementMap(
int o, map<LibUtilities::ShapeType, DerivUtilSharedPtr> derMap)
{
......
......@@ -172,6 +172,7 @@ void ProcessVarOpti::Process()
m_res = boost::shared_ptr<Residual>(new Residual);
m_res->val = 1.0;
m_mesh->MakeOrder(m_mesh->m_nummode - 1,
LibUtilities::eGaussLobattoLegendre);
......@@ -183,6 +184,7 @@ void ProcessVarOpti::Process()
map<LibUtilities::ShapeType, DerivUtilSharedPtr> derivUtils =
BuildDerivUtil(intOrder);
GetElementMap(intOrder, derivUtils);
m_res->startInv = 0;
......@@ -200,9 +202,10 @@ void ProcessVarOpti::Process()
elLock = GetLockedElements(m_config["region"].as<NekDouble>());
}
vector<vector<NodeSharedPtr> > freenodes = GetColouredNodes(elLock);
vector<vector<NodeOptiSharedPtr> > optiNodes;
// turn the free nodes into optimisable objects with all required data
set<int> check;
for (int i = 0; i < freenodes.size(); i++)
......
......@@ -104,8 +104,11 @@ private:
void GetElementMap(
int o, std::map<LibUtilities::ShapeType, DerivUtilSharedPtr> derMap);
std::vector<ElUtilSharedPtr> GetLockedElements(NekDouble thres);
std::vector<std::vector<NodeSharedPtr> > CreateColoursets(
std::vector<NodeSharedPtr> remain);
std::vector<std::vector<NodeSharedPtr> > GetColouredNodes(
std::vector<ElUtilSharedPtr> elLock);
void RemoveLinearCurvature();
NodeElMap m_nodeElMap;
......
......@@ -11,7 +11,7 @@
<regex>Worst at end: (-?\d+(?:\.\d*)?(?:[eE][+\-]?\d+)?)</regex>
<matches>
<match>
<field id="0" tolerance="5e-2">7.563667e-01</field>
<field id="0" tolerance="5e-2">7.042539e-01</field>
</match>
</matches>
</metric>
......
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