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

Merge branch 'master' into 'fix/mpi-cmake'

# Conflicts:
#   CHANGELOG.md
parents 3cf92905 271579a6
......@@ -34,8 +34,12 @@ v4.4.0
- Enabled MUMPS support in PETSc if a Fortran compiler was found and added 3D
support to the Helmholtz smoother used e.g. in FieldConverts C0Projection
module (!714)
- Fix bug in `Vmath::FillWhiteNoise` which caused `ForcingNoise` to have
a repeated pattern (!718)
- Fix bug in the calculation of the RHS magnitude in CG solver (!721)
- Fix bug in MPI detection for recent CMake on OS X (!725)
- Fix bug in CMake Homebrew and MacPorts detection for OS X (!729)
- Fix bug in FieldUtils when using half mode expansions (!734)
**ADRSolver:**
- Add a projection equation system for C^0 projections (!675)
......@@ -84,6 +88,11 @@ v4.4.0
(!712)
- 2D to 3D mesh extrusion module (!715)
- Add new two-dimensional mesher from NACA code or step file (!720)
- Fix inverted boundary layer in 2D (!736)
- More sensible element sizing with boundary layers in 2D (!736)
- 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)
**FieldConvert:**
- Move all modules to a new library, FieldUtils, to support post-processing
......
......@@ -190,7 +190,10 @@ state of the solution fields at at given timestep. This can subsequently be used
for restarting the simulation or examining time-dependent behaviour. This
produces a sequence of files, by default named \inltt{session\_*.chk}, where
\inltt{*} is replaced by a counter. The initial condition is written to
\inltt{session\_0.chk}.
\inltt{session\_0.chk}. Existing files are not overwritten, but renamed to e.g.
\inltt{session\_0.bak0.chk}. In case this file already exists, too, the \inltt{chk}-file
is renamed to \inltt{session\_0.bak*.chk} and so on.
\begin{notebox}
This functionality is equivalent to setting the \inltt{IO\_CheckSteps}
......
......@@ -281,7 +281,9 @@ struct Field
// Define Homogeneous expansion
int nplanes;
NekDouble lz;
LibUtilities::BasisType btype;
LibUtilities::BasisType btype;
LibUtilities::PointsType ptype =
LibUtilities::eFourierEvenlySpaced;
if (fldfilegiven)
{
......@@ -299,6 +301,11 @@ struct Field
nplanes = 4;
}
}
else if (btype == LibUtilities::eFourierHalfModeRe &&
nplanes == 1)
{
ptype = LibUtilities::ePolyEvenlySpaced;
}
}
else
{
......@@ -310,7 +317,7 @@ struct Field
// Choose points to be at evenly spaced points at
// nplanes points
const LibUtilities::PointsKey Pkey(
nplanes, LibUtilities::eFourierEvenlySpaced);
nplanes, ptype);
const LibUtilities::BasisKey Bkey(btype, nplanes, Pkey);
......
......@@ -100,19 +100,19 @@ class ThreadJob
{
public:
/// Base constructor
ThreadJob();
LIB_UTILITIES_EXPORT ThreadJob();
/// Base destructor.
virtual ~ThreadJob();
LIB_UTILITIES_EXPORT virtual ~ThreadJob();
/**
* This method will be called when the task is loaded
* onto a worker thread and is ready to run. When Run
* has finished this instance will be destructed.
*/
virtual void Run() = 0;
LIB_UTILITIES_EXPORT virtual void Run() = 0;
/// Set number of worker threads.
void SetWorkerNum(unsigned int num);
LIB_UTILITIES_EXPORT void SetWorkerNum(unsigned int num);
protected:
/**
......@@ -168,7 +168,7 @@ class ThreadManager : public boost::enable_shared_from_this<ThreadManager>
{
public:
/// Destructor.
virtual ~ThreadManager();
LIB_UTILITIES_EXPORT virtual ~ThreadManager();
/**
* @brief Pass a list of tasklets to the master queue.
* @param joblist Vector of ThreadJob pointers.
......@@ -181,7 +181,7 @@ class ThreadManager : public boost::enable_shared_from_this<ThreadManager>
*
* @see SchedType
*/
virtual void QueueJobs(std::vector<ThreadJob*>& joblist) = 0;
LIB_UTILITIES_EXPORT virtual void QueueJobs(std::vector<ThreadJob*>& joblist) = 0;
/**
* @brief Pass a single job to the master queue.
* @param job A pointer to a ThreadJob subclass.
......@@ -190,14 +190,14 @@ class ThreadManager : public boost::enable_shared_from_this<ThreadManager>
* issue then suspend the workers with SetNumWorkers(0) until the jobs
* are queued.
*/
virtual void QueueJob(ThreadJob* job) = 0;
LIB_UTILITIES_EXPORT virtual void QueueJob(ThreadJob* job) = 0;
/**
* @brief Return the number of active workers.
*
* Active workers are threads that are either running jobs
* or are waiting for jobs to be queued.
*/
virtual unsigned int GetNumWorkers() = 0;
LIB_UTILITIES_EXPORT virtual unsigned int GetNumWorkers() = 0;
/**
* @brief Returns the worker number of the executing thread.
*
......@@ -216,7 +216,7 @@ class ThreadManager : public boost::enable_shared_from_this<ThreadManager>
*
* Returns 0 if called by non-thread.
*/
virtual unsigned int GetWorkerNum() = 0;
LIB_UTILITIES_EXPORT virtual unsigned int GetWorkerNum() = 0;
/**
* @brief Sets the number of active workers.
* @param num The number of active workers.
......@@ -227,18 +227,18 @@ class ThreadManager : public boost::enable_shared_from_this<ThreadManager>
* If num is greater than the maximum allowed number of active workers,
* then the maximum value will be used instead.
*/
virtual void SetNumWorkers(const unsigned int num) = 0;
LIB_UTILITIES_EXPORT virtual void SetNumWorkers(const unsigned int num) = 0;
/**
* @brief Sets the number of active workers to the maximum.
*
* Sets the number of active workers to the maximum available.
*/
virtual void SetNumWorkers() = 0;
LIB_UTILITIES_EXPORT virtual void SetNumWorkers() = 0;
/**
* @brief Gets the maximum available number of threads.
* @return The maximum number of workers.
*/
virtual unsigned int GetMaxNumWorkers() = 0;
LIB_UTILITIES_EXPORT virtual unsigned int GetMaxNumWorkers() = 0;
/**
* @brief Waits until all queued jobs are finished.
*
......@@ -261,7 +261,7 @@ class ThreadManager : public boost::enable_shared_from_this<ThreadManager>
* number of worker threads, implementations should increase the number
* of active workers by 1 on entering Wait().
*/
virtual void Wait() = 0;
LIB_UTILITIES_EXPORT virtual void Wait() = 0;
/**
* @brief Controls how many jobs are sent to each worker at a time.
*
......@@ -271,17 +271,17 @@ class ThreadManager : public boost::enable_shared_from_this<ThreadManager>
* @see SchedType
* @see SetSchedType()
*/
virtual void SetChunkSize(unsigned int chnk) = 0;
LIB_UTILITIES_EXPORT virtual void SetChunkSize(unsigned int chnk) = 0;
/**
* @brief Sets the current scheduling algorithm.
* @see SetChunkSize()
*/
virtual void SetSchedType(SchedType s) = 0;
LIB_UTILITIES_EXPORT virtual void SetSchedType(SchedType s) = 0;
/**
* @brief Indicates whether the code is in a worker thread or not.
* @return True if the caller is in a worker thread.
*/
virtual bool InThread() = 0;
LIB_UTILITIES_EXPORT virtual bool InThread() = 0;
/**
* @brief A calling threads holds until all active threads call this
* method.
......@@ -294,16 +294,16 @@ class ThreadManager : public boost::enable_shared_from_this<ThreadManager>
* is altered after a thread has called this method. It is only safe to
* call SetNumWorkers() when no threads are holding.
*/
virtual void Hold() = 0;
LIB_UTILITIES_EXPORT virtual void Hold() = 0;
/**
* @brief Returns a description of the type of threading.
*
* E.g. "Threading with Boost"
*/
virtual const std::string& GetType() const = 0;
LIB_UTILITIES_EXPORT virtual const std::string& GetType() const = 0;
/// ThreadManager implementation.
virtual bool IsInitialised();
LIB_UTILITIES_EXPORT virtual bool IsInitialised();
inline int GetThrFromPartition(int pPartition)
{
......@@ -354,17 +354,17 @@ class ThreadMaster
THREADMANAGER_MAX
};
/// Constructor
ThreadMaster();
LIB_UTILITIES_EXPORT ThreadMaster();
/// Destructor
~ThreadMaster();
LIB_UTILITIES_EXPORT ~ThreadMaster();
/// Sets what ThreadManagers will be created in CreateInstance.
void SetThreadingType(const std::string &p_type);
LIB_UTILITIES_EXPORT void SetThreadingType(const std::string &p_type);
/// Gets the ThreadManager associated with string s.
ThreadManagerSharedPtr& GetInstance(const ThreadManagerName t);
LIB_UTILITIES_EXPORT ThreadManagerSharedPtr& GetInstance(const ThreadManagerName t);
/// Creates an instance of a ThreadManager (which one is determined by
/// a previous call to SetThreadingType) and associates it with
/// the string s.
ThreadManagerSharedPtr CreateInstance(const ThreadManagerName t,
LIB_UTILITIES_EXPORT ThreadManagerSharedPtr CreateInstance(const ThreadManagerName t,
unsigned int nThr);
};
......
......@@ -135,30 +135,44 @@ namespace Vmath
template LIB_UTILITIES_EXPORT Nektar::NekDouble ran2 (long* idum);
/// \brief Fills a vector with white noise.
template<class T> void FillWhiteNoise( int n, const T eps, T *x, const int incx, int outseed)
template<class T> void FillWhiteNoise( int n, const T eps, T *x,
const int incx, int outseed)
{
// Protect the static vars here and in ran2
boost::mutex::scoped_lock l(mutex);
// Define static variables for generating random numbers
static int iset = 0;
static T gset;
static long seed = 0;
// Bypass seed if outseed was specified
if( outseed != 9999)
{
seed = long(outseed);
}
while( n-- )
{
static int iset = 0;
static T gset;
long seed = long(outseed);
T fac, rsq, v1, v2;
if (iset == 0) {
do {
v1 = 2.0 * ran2<T> (&seed) - 1.0;
v2 = 2.0 * ran2<T> (&seed) - 1.0;
rsq = v1*v1 + v2*v2;
} while (rsq >= 1.0 || rsq == 0.0);
fac = sqrt(-2.0 * log (rsq) / rsq);
gset = v1 * fac;
iset = 1;
*x = eps * v2 * fac;
} else {
iset = 0;
*x = eps * gset;
if (iset == 0)
{
do
{
v1 = 2.0 * ran2<T> (&seed) - 1.0;
v2 = 2.0 * ran2<T> (&seed) - 1.0;
rsq = v1*v1 + v2*v2;
} while (rsq >= 1.0 || rsq == 0.0);
fac = sqrt(-2.0 * log (rsq) / rsq);
gset = v1 * fac;
iset = 1;
*x = eps * v2 * fac;
}
else
{
iset = 0;
*x = eps * gset;
}
x += incx;
}
......
......@@ -54,7 +54,8 @@ namespace Vmath
/// \brief Fills a vector with white noise.
template<class T> LIB_UTILITIES_EXPORT void FillWhiteNoise( int n, const T eps, T *x, const int incx, int seed = 0);
template<class T> LIB_UTILITIES_EXPORT void FillWhiteNoise(
int n, const T eps, T *x, const int incx, int seed = 9999);
/// \brief Multiply vector z = x*y
template<class T> LIB_UTILITIES_EXPORT void Vmul( int n, const T *x, const int incx, const T *y,
......
......@@ -53,7 +53,9 @@
Fill(n,alpha,&x[0],incx);
}
template<class T> void FillWhiteNoise( int n, const T eps, Array<OneD, T> &x, const int incx, int outseed)
template<class T> void FillWhiteNoise(
int n, const T eps, Array<OneD, T> &x,
const int incx, int outseed = 9999)
{
ASSERTL1(n*incx <= x.num_elements()+x.GetOffset(),"Out of bounds");
......
......@@ -418,6 +418,8 @@ namespace Nektar
Set_Rhs_Magnitude(inGlob);
}
m_totalIterations = 0;
// If input residual is less than tolerance skip solve.
if (eps < m_tolerance * m_tolerance * m_rhs_magnitude)
{
......@@ -432,7 +434,6 @@ namespace Nektar
return;
}
m_totalIterations = 1;
m_precon->DoPreconditioner(r_A, tmp = w_A + nDir);
v_DoMatrixMultiply(w_A, s_A);
......@@ -451,11 +452,12 @@ namespace Nektar
vComm->AllReduce(vExchange, Nektar::LibUtilities::ReduceSum);
rho = vExchange[0];
mu = vExchange[1];
min_resid = m_rhs_magnitude;
beta = 0.0;
alpha = rho/mu;
rho = vExchange[0];
mu = vExchange[1];
min_resid = m_rhs_magnitude;
beta = 0.0;
alpha = rho/mu;
m_totalIterations = 1;
// Continue until convergence
while (true)
......@@ -541,8 +543,12 @@ namespace Nektar
void GlobalLinSysIterative::Set_Rhs_Magnitude(
const NekVector<NekDouble> &pIn)
{
Array<OneD, NekDouble> vExchange(1);
vExchange[0] = Vmath::Dot(pIn.GetDimension(),&pIn[0],&pIn[0]);
Array<OneD, NekDouble> vExchange(1, 0.0);
if (m_map.num_elements() > 0)
{
vExchange[0] = Vmath::Dot2(pIn.GetDimension(),
&pIn[0],&pIn[0],&m_map[0]);
}
m_expList.lock()->GetComm()->GetRowComm()->AllReduce(
vExchange, Nektar::LibUtilities::ReduceSum);
......
......@@ -54,7 +54,7 @@ Generator2D::Generator2D(MeshSharedPtr m) : ProcessModule(m)
m_config["blcurves"] =
ConfigOption(false, "0", "Generate parallelograms on these curves");
m_config["blthick"] =
ConfigOption(false, "0", "Parallelogram layer thickness");
ConfigOption(false, "0.0", "Parallelogram layer thickness");
}
Generator2D::~Generator2D()
......@@ -71,6 +71,9 @@ void Generator2D::Process()
m_mesh->m_numNodes = m_mesh->m_cad->GetNumVerts();
m_thickness_ID =
m_thickness.DefineFunction("x y z", m_config["blthick"].as<string>());
// linear mesh all curves
for (int i = 1; i <= m_mesh->m_cad->GetNumCurve(); i++)
{
......@@ -80,12 +83,25 @@ void Generator2D::Process()
"Curve progress");
}
m_curvemeshes[i] =
MemoryManager<CurveMesh>::AllocateSharedPtr(i, m_mesh);
vector<unsigned int>::iterator f =
find(m_blCurves.begin(), m_blCurves.end(), i);
if (f == m_blCurves.end())
{
m_curvemeshes[i] =
MemoryManager<CurveMesh>::AllocateSharedPtr(i, m_mesh);
}
else
{
m_curvemeshes[i] = MemoryManager<CurveMesh>::AllocateSharedPtr(
i, m_mesh, m_config["blthick"].as<string>());
}
m_curvemeshes[i]->Mesh();
}
////////////////////////////////////////
EdgeSet::iterator it;
for (it = m_mesh->m_edgeSet.begin(); it != m_mesh->m_edgeSet.end(); it++)
{
......@@ -108,6 +124,14 @@ void Generator2D::Process()
m_mesh->m_element[1].push_back(E2);
}
for (int i = 1; i <= m_mesh->m_cad->GetNumSurf(); i++)
{
m_facemeshes[i] = MemoryManager<FaceMesh>::AllocateSharedPtr(
i, m_mesh, m_curvemeshes, 100);
m_facemeshes[i]->OrientateCurves();
}
if (m_config["blcurves"].beenSet)
{
// we need to do the boundary layer generation in a face by face basis
......@@ -120,15 +144,11 @@ void Generator2D::Process()
// orienation
for (int i = 1; i <= m_mesh->m_cad->GetNumSurf(); i++)
{
m_facemeshes[i] = MemoryManager<FaceMesh>::AllocateSharedPtr(
i, m_mesh, m_curvemeshes, 100);
m_facemeshes[i]->OrientateCurves();
MakeBL(i, m_facemeshes[i]->GetEdges());
}
}
//m_mesh->m_element[1].clear();
// m_mesh->m_element[1].clear();
if (m_mesh->m_verbose)
{
......@@ -136,7 +156,9 @@ void Generator2D::Process()
}
// linear mesh all surfaces
for (int i = 1; i <= m_mesh->m_cad->GetNumSurf(); i++)
map<int, FaceMeshSharedPtr>::iterator fit;
int i = 1;
for (fit = m_facemeshes.begin(); fit != m_facemeshes.end(); fit++)
{
if (m_mesh->m_verbose)
{
......@@ -144,10 +166,16 @@ void Generator2D::Process()
"Face progress");
}
m_facemeshes[i] = MemoryManager<FaceMesh>::AllocateSharedPtr(
i, m_mesh, m_curvemeshes, 100);
if (m_config["blcurves"].beenSet)
{
// for bl surfaces orientate curves needs to be run again to
// push nodes to the edges of the system
fit->second->ResetCurvemeshes(m_curvemeshes);
fit->second->OrientateCurves();
}
m_facemeshes[i]->Mesh();
fit->second->Mesh();
i++;
}
ProcessVertices();
......@@ -169,7 +197,6 @@ void Generator2D::MakeBLPrep()
// identify the nodes which will become the boundary layer.
ParseUtils::GenerateSeqVector(m_config["blcurves"].as<string>().c_str(),
m_blCurves);
m_thickness = m_config["blthick"].as<NekDouble>();
for (vector<unsigned>::iterator it = m_blCurves.begin();
it != m_blCurves.end(); ++it)
......@@ -183,14 +210,14 @@ void Generator2D::MakeBLPrep()
}
}
void Generator2D::MakeBL(int faceid, vector<EdgeLoop> e)
void Generator2D::MakeBL(int faceid, vector<EdgeLoopSharedPtr> e)
{
map<int, int> edgeToOrient;
for (vector<EdgeLoop>::iterator lit = e.begin(); lit != e.end(); ++lit)
for (int j = 0; j < e.size(); j++)
{
for (int i = 0; i < lit->edges.size(); ++i)
for (int i = 0; i < e[j]->edges.size(); ++i)
{
edgeToOrient[lit->edges[i]->GetId()] = lit->edgeo[i];
edgeToOrient[e[j]->edges[i]->GetId()] = e[j]->edgeo[i];
}
}
......@@ -199,7 +226,7 @@ void Generator2D::MakeBL(int faceid, vector<EdgeLoop> e)
int eid = 0;
for (vector<unsigned>::iterator it = m_blCurves.begin();
it != m_blCurves.end(); ++it)
it != m_blCurves.end(); ++it)
{
int edgeo = edgeToOrient[*it];
......@@ -207,17 +234,17 @@ void Generator2D::MakeBL(int faceid, vector<EdgeLoop> e)
// on each !!!EDGE!!! calculate a normal
// always to the left unless edgeo is 1
for(int j = 0; j < es.size(); j++)
for (int j = 0; j < es.size(); j++)
{
es[j]->m_id = eid++;
Array<OneD, NekDouble> p1 = (edgeo == 0) ? es[j]->m_n1->GetLoc()
: es[j]->m_n2->GetLoc();
Array<OneD, NekDouble> p2 = (edgeo == 0) ? es[j]->m_n2->GetLoc()
: es[j]->m_n1->GetLoc();
Array<OneD, NekDouble> p1 =
(edgeo == 0) ? es[j]->m_n1->GetLoc() : es[j]->m_n2->GetLoc();
Array<OneD, NekDouble> p2 =
(edgeo == 0) ? es[j]->m_n2->GetLoc() : es[j]->m_n1->GetLoc();
Array<OneD, NekDouble> n(2);
n[0] = p1[1] - p2[1];
n[1] = p2[0] - p1[0];
NekDouble mag = sqrt(n[0]*n[0]+n[1]*n[1]);
n[0] = p1[1] - p2[1];
n[1] = p2[0] - p1[0];
NekDouble mag = sqrt(n[0] * n[0] + n[1] * n[1]);
n[0] /= mag;
n[1] /= mag;
edgeNormals[es[j]->m_id] = n;
......@@ -226,7 +253,7 @@ void Generator2D::MakeBL(int faceid, vector<EdgeLoop> e)
map<NodeSharedPtr, NodeSharedPtr> nodeNormals;
map<NodeSharedPtr, vector<EdgeSharedPtr> >::iterator it;
for(it = m_nodesToEdge.begin(); it != m_nodesToEdge.end(); it++)
for (it = m_nodesToEdge.begin(); it != m_nodesToEdge.end(); it++)
{
Array<OneD, NekDouble> n(3);
ASSERTL0(it->second.size() == 2,
......@@ -234,39 +261,42 @@ void Generator2D::MakeBL(int faceid, vector<EdgeLoop> e)
Array<OneD, NekDouble> n1 = edgeNormals[it->second[0]->m_id];
Array<OneD, NekDouble> n2 = edgeNormals[it->second[1]->m_id];
n[0] = (n1[0] + n2[0]) / 2.0;
n[1] = (n1[1] + n2[1]) / 2.0;
NekDouble mag = sqrt(n[0]*n[0]+n[1]*n[1]);
n[0] = (n1[0] + n2[0]) / 2.0;
n[1] = (n1[1] + n2[1]) / 2.0;
NekDouble mag = sqrt(n[0] * n[0] + n[1] * n[1]);
n[0] /= mag;
n[1] /= mag;
n[0] = n[0] * m_thickness + it->first->m_x;
n[1] = n[1] * m_thickness + it->first->m_y;
NekDouble t = m_thickness.Evaluate(m_thickness_ID, it->first->m_x,
it->first->m_y, 0.0, 0.0);
n[0] = n[0] * t + it->first->m_x;
n[1] = n[1] * t + it->first->m_y;
n[2] = 0.0;
NodeSharedPtr nn = boost::shared_ptr<Node>(
new Node(m_mesh->m_numNodes++, n[0], n[1], 0.0));
CADSurfSharedPtr s = m_mesh->m_cad->GetSurf(faceid);
Array<OneD, NekDouble> uv = s->locuv(n);
nn->SetCADSurf(faceid,s,uv);
nn->SetCADSurf(faceid, s, uv);
nodeNormals[it->first] = nn;
}
for (vector<unsigned>::iterator it = m_blCurves.begin();
it != m_blCurves.end(); ++it)
it != m_blCurves.end(); ++it)
{
int edgeo = edgeToOrient[*it];
vector<NodeSharedPtr> ns = m_curvemeshes[*it]->GetMeshPoints();
vector<NodeSharedPtr> newNs;
for(int i = 0; i < ns.size(); i++)
for (int i = 0; i < ns.size(); i++)
{
newNs.push_back(nodeNormals[ns[i]]);
}
m_curvemeshes[*it] = MemoryManager<CurveMesh>::AllocateSharedPtr(
*it, m_mesh, newNs);
m_curvemeshes[*it] =
MemoryManager<CurveMesh>::AllocateSharedPtr(*it, m_mesh, newNs);