Nektar++
Public Member Functions | Static Public Member Functions | Static Public Attributes | Protected Member Functions | Protected Attributes | Friends | List of all members
Nektar::SolverUtils::FilterHistoryPoints Class Reference

#include <FilterHistoryPoints.h>

Inheritance diagram for Nektar::SolverUtils::FilterHistoryPoints:
[legend]

Public Member Functions

SOLVER_UTILS_EXPORT FilterHistoryPoints (const LibUtilities::SessionReaderSharedPtr &pSession, const std::shared_ptr< EquationSystem > &pEquation, const ParamMap &pParams)
 
SOLVER_UTILS_EXPORT ~FilterHistoryPoints () override
 
- Public Member Functions inherited from Nektar::SolverUtils::Filter
SOLVER_UTILS_EXPORT Filter (const LibUtilities::SessionReaderSharedPtr &pSession, const std::shared_ptr< EquationSystem > &pEquation)
 
virtual SOLVER_UTILS_EXPORT ~Filter ()
 
SOLVER_UTILS_EXPORT void Initialise (const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
 
SOLVER_UTILS_EXPORT void Update (const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
 
SOLVER_UTILS_EXPORT void Finalise (const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
 
SOLVER_UTILS_EXPORT bool IsTimeDependent ()
 
SOLVER_UTILS_EXPORT std::string SetupOutput (const std::string ext, const ParamMap &pParams)
 
SOLVER_UTILS_EXPORT std::string SetupOutput (const std::string ext, const std::string inname)
 
SOLVER_UTILS_EXPORT void SetUpdateOnInitialise (bool flag)
 

Static Public Member Functions

static FilterSharedPtr create (const LibUtilities::SessionReaderSharedPtr &pSession, const std::shared_ptr< EquationSystem > &pEquation, const std::map< std::string, std::string > &pParams)
 Creates an instance of this class. More...
 

Static Public Attributes

static std::string className
 Name of the class. More...
 

Protected Member Functions

SOLVER_UTILS_EXPORT void v_Initialise (const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time) override
 
SOLVER_UTILS_EXPORT void v_Update (const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time) override
 
SOLVER_UTILS_EXPORT void v_Finalise (const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time) override
 
SOLVER_UTILS_EXPORT bool v_IsTimeDependent () override
 
bool GetPoint (Array< OneD, NekDouble > gloCoord, int I)
 
SOLVER_UTILS_EXPORT void v_WriteData (const int &rank, const Array< OneD, NekDouble > &data, const int &numFields, const NekDouble &time)
 
- Protected Member Functions inherited from Nektar::SolverUtils::Filter
virtual void v_Initialise (const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)=0
 
virtual void v_Update (const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)=0
 
virtual void v_Finalise (const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)=0
 
virtual bool v_IsTimeDependent ()=0
 
virtual SOLVER_UTILS_EXPORT std::string v_SetupOutput (const std::string ext, const ParamMap &pParams)
 
virtual SOLVER_UTILS_EXPORT std::string v_SetupOutput (const std::string ext, const std::string inname)
 

Protected Attributes

Array< OneD, Array< OneD, const NekDouble > > m_historyPoints
 
size_t m_historyPointsSize = 0
 
unsigned int m_index = 0
 
unsigned int m_outputFrequency
 
int m_outputPlane
 plane to take history point from if using a homogeneous1D expansion More...
 
std::vector< int > m_planeIDs
 
bool m_isHomogeneous1D
 
bool m_waveSpace
 
std::string m_outputFile
 
std::ofstream m_outputStream
 
std::stringstream m_historyPointStream
 
std::list< std::tuple< Array< OneD, const NekDouble >, Array< OneD, const NekDouble >, int, int > > m_historyList
 
std::map< LibUtilities::PtsType, Array< OneD, NekDouble > > m_pointDatMap
 
std::map< LibUtilities::PtsType, Array< OneD, int > > m_pointNumMap
 
unsigned int m_outputIndex = 0
 
bool m_outputOneFile
 
bool m_adaptive
 
- Protected Attributes inherited from Nektar::SolverUtils::Filter
LibUtilities::SessionReaderSharedPtr m_session
 
const std::weak_ptr< EquationSystemm_equ
 
bool m_updateOnInitialise = true
 

Friends

class MemoryManager< FilterHistoryPoints >
 

Additional Inherited Members

- Public Types inherited from Nektar::SolverUtils::Filter
typedef std::map< std::string, std::string > ParamMap
 

Detailed Description

Definition at line 45 of file FilterHistoryPoints.h.

Constructor & Destructor Documentation

◆ FilterHistoryPoints()

Nektar::SolverUtils::FilterHistoryPoints::FilterHistoryPoints ( const LibUtilities::SessionReaderSharedPtr pSession,
const std::shared_ptr< EquationSystem > &  pEquation,
const ParamMap pParams 
)

Definition at line 53 of file FilterHistoryPoints.cpp.

56 : Filter(pSession, pEquation)
57{
58 // OutputFile
59 std::string ext = ".his";
60 m_outputFile = Filter::SetupOutput(ext, pParams);
61
62 // OutputFrequency
63 auto it = pParams.find("OutputFrequency");
64 if (it == pParams.end())
65 {
67 }
68 else
69 {
70 LibUtilities::Equation equ(m_session->GetInterpreter(), it->second);
71 m_outputFrequency = round(equ.Evaluate());
72 }
73
74 // output all data into one file
75 it = pParams.find("OutputOneFile");
76 if (it == pParams.end())
77 {
78 m_outputOneFile = true;
79 }
80 else
81 {
82 std::string sOption = it->second.c_str();
83 m_outputOneFile = (boost::iequals(sOption, "true")) ||
84 (boost::iequals(sOption, "yes"));
85 }
86
87 // OutputPlane
88 m_session->MatchSolverInfo("Homogeneous", "1D", m_isHomogeneous1D, false);
90 {
91 it = pParams.find("OutputPlane");
92 if (it == pParams.end())
93 {
94 m_outputPlane = -1;
95 }
96 else
97 {
98 LibUtilities::Equation equ(m_session->GetInterpreter(), it->second);
99 m_outputPlane = round(equ.Evaluate());
100 }
101
102 it = pParams.find("WaveSpace");
103 if (it == pParams.end())
104 {
105 m_waveSpace = false;
106 }
107 else
108 {
109 std::string sOption = it->second.c_str();
110 m_waveSpace = (boost::iequals(sOption, "true")) ||
111 (boost::iequals(sOption, "yes"));
112 }
113 }
114
115 // Points
116 if (pParams.end() != (it = pParams.find("Points")))
117 {
118 m_pointNumMap[LibUtilities::ePtsFile] = Array<OneD, int>(1, 0);
119 m_historyPointStream.str(it->second);
120 }
121 else if (pParams.end() != (it = pParams.find("line")))
122 {
123 vector<NekDouble> values;
124 ASSERTL0(ParseUtils::GenerateVector(it->second, values),
125 "Failed to interpret line string");
126
127 ASSERTL0(values.size() > 2, "line string should contain 2*Dim+1 values "
128 "N,x0,y0,z0,x1,y1,z1");
129
130 double tmp;
131 ASSERTL0(std::modf(values[0], &tmp) == 0.0, "N is not an integer");
132 ASSERTL0(values[0] > 1, "N is not a valid number");
133
134 int dim = (values.size() - 1) / 2;
135 int npts = values[0];
136
137 Array<OneD, int> num(1, npts);
138 Array<OneD, NekDouble> delta(6, 0.);
139 for (int i = 0; i < dim; ++i)
140 {
141 delta[i + 0] = values[i + 1];
142 delta[i + 3] = (values[dim + i + 1] - values[i + 1]) / (npts - 1);
143 }
146 }
147 else if (pParams.end() != (it = pParams.find("plane")))
148 {
149 vector<NekDouble> values;
150 ASSERTL0(ParseUtils::GenerateVector(it->second, values),
151 "Failed to interpret plane string");
152
153 ASSERTL0(values.size() > 9,
154 "plane string should contain 4 Dim+2 values "
155 "N1,N2,x0,y0,z0,x1,y1,z1,x2,y2,z2,x3,y3,z3");
156
157 double tmp;
158 ASSERTL0(std::modf(values[0], &tmp) == 0.0, "N1 is not an integer");
159 ASSERTL0(std::modf(values[1], &tmp) == 0.0, "N2 is not an integer");
160
161 ASSERTL0(values[0] > 1, "N1 is not a valid number");
162 ASSERTL0(values[1] > 1, "N2 is not a valid number");
163
164 int dim = (values.size() - 2) / 4;
165
166 Array<OneD, int> npts(3);
167 npts[0] = values[0];
168 npts[1] = values[0] * values[1];
169 npts[2] = values[1];
170
171 Array<OneD, NekDouble> delta(12, 0.);
172 for (int i = 0; i < dim; ++i)
173 {
174 delta[i + 0] = values[2 + i];
175 delta[i + 3] = values[2 + 3 * dim + i];
176 delta[i + 6] = (values[2 + 1 * dim + i] - values[2 + 0 * dim + i]) /
177 (values[0] - 1);
178 delta[i + 9] = (values[2 + 2 * dim + i] - values[2 + 3 * dim + i]) /
179 (values[0] - 1);
180 }
183 }
184 else if (pParams.end() != (it = pParams.find("box")))
185 {
186 vector<NekDouble> values;
187 ASSERTL0(ParseUtils::GenerateVector(it->second, values),
188 "Failed to interpret box string");
189
190 ASSERTL0(values.size() == 9, "box string should contain 9 values "
191 "N1,N2,N3,xmin,xmax,ymin,ymax,zmin,zmax");
192
193 int dim = 3;
194 Array<OneD, int> npts(3);
195 npts[0] = values[0];
196 npts[1] = values[0] * values[1];
197 npts[2] = values[0] * values[1] * values[2];
198
199 Array<OneD, NekDouble> delta(6, 0.);
200 for (int i = 0; i < dim; ++i)
201 {
202 delta[i + 0] = values[3 + 2 * i];
203 delta[i + 3] =
204 (values[4 + 2 * i] - values[3 + 2 * i]) / (values[i] - 1);
205 }
208 }
209 else
210 {
211 ASSERTL0(false, "Missing parameter 'Points'.");
212 }
213}
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
static bool GenerateVector(const std::string &str, std::vector< T > &out)
Takes a comma-separated string and converts it to entries in a vector.
Definition: ParseUtils.cpp:130
std::map< LibUtilities::PtsType, Array< OneD, int > > m_pointNumMap
int m_outputPlane
plane to take history point from if using a homogeneous1D expansion
std::map< LibUtilities::PtsType, Array< OneD, NekDouble > > m_pointDatMap
SOLVER_UTILS_EXPORT std::string SetupOutput(const std::string ext, const ParamMap &pParams)
Definition: Filter.h:139
LibUtilities::SessionReaderSharedPtr m_session
Definition: Filter.h:93
SOLVER_UTILS_EXPORT Filter(const LibUtilities::SessionReaderSharedPtr &pSession, const std::shared_ptr< EquationSystem > &pEquation)

References ASSERTL0, Nektar::LibUtilities::ePtsBox, Nektar::LibUtilities::ePtsFile, Nektar::LibUtilities::ePtsLine, Nektar::LibUtilities::ePtsPlane, Nektar::LibUtilities::Equation::Evaluate(), Nektar::ParseUtils::GenerateVector(), m_historyPointStream, m_isHomogeneous1D, m_outputFile, m_outputFrequency, m_outputOneFile, m_outputPlane, m_pointDatMap, m_pointNumMap, Nektar::SolverUtils::Filter::m_session, m_waveSpace, and Nektar::SolverUtils::Filter::SetupOutput().

◆ ~FilterHistoryPoints()

Nektar::SolverUtils::FilterHistoryPoints::~FilterHistoryPoints ( )
override

Definition at line 218 of file FilterHistoryPoints.cpp.

219{
220}

Member Function Documentation

◆ create()

static FilterSharedPtr Nektar::SolverUtils::FilterHistoryPoints::create ( const LibUtilities::SessionReaderSharedPtr pSession,
const std::shared_ptr< EquationSystem > &  pEquation,
const std::map< std::string, std::string > &  pParams 
)
inlinestatic

Creates an instance of this class.

Definition at line 51 of file FilterHistoryPoints.h.

55 {
58 pSession, pEquation, pParams);
59 return p;
60 }
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
std::shared_ptr< Filter > FilterSharedPtr
A shared pointer to a Driver object.
Definition: Filter.h:52

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), and CellMLToNektar.cellml_metadata::p.

◆ GetPoint()

bool Nektar::SolverUtils::FilterHistoryPoints::GetPoint ( Array< OneD, NekDouble gloCoord,
int  I 
)
protected

Definition at line 222 of file FilterHistoryPoints.cpp.

223{
225 {
226 m_historyPointStream >> gloCoord[0] >> gloCoord[1] >> gloCoord[2];
227 return !m_historyPointStream.fail();
228 }
230 {
232 {
233 return false;
234 }
235 Array<OneD, NekDouble> values = m_pointDatMap[LibUtilities::ePtsLine];
236 Array<OneD, NekDouble> delta =
238 for (int n = 0; n < 3; ++n)
239 {
240 gloCoord[n] = values[n] + I * delta[n];
241 }
242 return true;
243 }
245 {
247 {
248 return false;
249 }
250 Array<OneD, NekDouble> values = m_pointDatMap[LibUtilities::ePtsPlane];
251 Array<OneD, NekDouble> delta =
253 Array<OneD, int> i(2);
255 i[0] = I - i[1] * m_pointNumMap[LibUtilities::ePtsPlane][0];
256 double n1 = -1. + (NekDouble)m_pointNumMap[LibUtilities::ePtsPlane][2];
257 for (int n = 0; n < 3; ++n)
258 {
259 gloCoord[n] = (values[n] + i[0] * delta[n]) * (1.0 - i[1] / n1) +
260 (values[3 + n] + i[0] * delta[3 + n]) * (i[1] / n1);
261 }
262 return true;
263 }
264 else if (m_pointNumMap.count(LibUtilities::ePtsBox))
265 {
267 {
268 return false;
269 }
270 Array<OneD, NekDouble> values = m_pointDatMap[LibUtilities::ePtsBox];
271 Array<OneD, NekDouble> delta = m_pointDatMap[LibUtilities::ePtsBox] + 3;
272 Array<OneD, int> i(3);
274 I -= i[2] * m_pointNumMap[LibUtilities::ePtsBox][1];
276 i[0] = I - i[1] * m_pointNumMap[LibUtilities::ePtsBox][0];
277 for (int n = 0; n < 3; ++n)
278 {
279 gloCoord[n] = values[n] + i[n] * delta[n];
280 }
281 return true;
282 }
283 return false;
284}
double NekDouble

References Nektar::LibUtilities::ePtsBox, Nektar::LibUtilities::ePtsFile, Nektar::LibUtilities::ePtsLine, Nektar::LibUtilities::ePtsPlane, m_historyPointStream, m_pointDatMap, and m_pointNumMap.

Referenced by v_Initialise().

◆ v_Finalise()

void Nektar::SolverUtils::FilterHistoryPoints::v_Finalise ( const Array< OneD, const MultiRegions::ExpListSharedPtr > &  pFields,
const NekDouble time 
)
overrideprotectedvirtual

Implements Nektar::SolverUtils::Filter.

Definition at line 681 of file FilterHistoryPoints.cpp.

684{
685 if (pFields[0]->GetComm()->GetRank() == 0 && m_outputOneFile)
686 {
687 if (m_outputStream.is_open())
688 {
689 m_outputStream.close();
690 }
691 }
692}

References m_outputOneFile, and m_outputStream.

◆ v_Initialise()

void Nektar::SolverUtils::FilterHistoryPoints::v_Initialise ( const Array< OneD, const MultiRegions::ExpListSharedPtr > &  pFields,
const NekDouble time 
)
overrideprotectedvirtual

Implements Nektar::SolverUtils::Filter.

Definition at line 289 of file FilterHistoryPoints.cpp.

292{
293 ASSERTL0(!m_historyPointStream.fail(), "No history points in stream.");
294
295 LibUtilities::CommSharedPtr vComm = pFields[0]->GetComm();
296
297 // All processes read all the history points
298 Array<OneD, NekDouble> gloCoord(3, 0.0);
299 std::vector<Array<OneD, const NekDouble>> GloCoords;
300 const int coordim = pFields[0]->GetGraph()->GetSpaceDimension();
301
302 int i = -1;
303 while (GetPoint(gloCoord, ++i))
304 {
305 // Overwrite gloCoord[2] for 3DH1D using m_outputPlane
306 // if it is defined
308 {
309 int nplanes = pFields[0]->GetHomogeneousBasis()->GetZ().size();
310 NekDouble lhom = pFields[0]->GetHomoLen();
311 int plane;
312 if (m_outputPlane != -1)
313 {
314 plane = m_outputPlane;
315 }
316 else
317 {
318 // Pick plane near the point
319 plane = round((gloCoord[coordim] * nplanes) / lhom);
320 }
321 if (m_waveSpace)
322 {
323 m_planeIDs.push_back(plane);
324 }
325 NekDouble Z = (pFields[0]->GetHomogeneousBasis()->GetZ())[plane];
326 Z = (Z + 1) * lhom / 2;
327 if (fabs(gloCoord[coordim] - Z) >
329 vComm->GetRank() == 0)
330 {
331 cout << "Resetting History point from z = " << gloCoord[coordim]
332 << " to z = " << Z << endl;
333 }
334 gloCoord[coordim] = Z;
335 }
336
337 // Save a copy of the global coordinate
338 GloCoords.push_back(
339 Array<OneD, const NekDouble>(gloCoord.size(), gloCoord.data()));
340 }
341
342 // Check which history points this process is responsible for
343 const int vRank = vComm->GetRank();
344 const int vColumnRank = vComm->GetColumnComm()->GetRank();
345 const int vRowRank = vComm->GetRowComm()->GetRank();
346 m_historyPointsSize = GloCoords.size();
347 Array<OneD, int> procList(m_historyPointsSize, -1);
348 Array<OneD, int> idList(m_historyPointsSize, -1);
349
350 Array<OneD, NekDouble> locCoord(pFields[0]->GetShapeDimension(), 0.);
351 std::vector<Array<OneD, const NekDouble>> LocCoords;
352
353 // For each history point, find the element that contains it. If no element
354 // is available on this process, idList[i] = -1 and procList[i] = -1
355 for (i = 0; i < m_historyPointsSize; ++i)
356 {
357 // Create a copy of the global coordinate
358 gloCoord = GloCoords[i];
359
360 // Determine the expansion and local coordinates
362 {
363 idList[i] = pFields[0]->GetPlane(0)->GetExpIndex(
364 gloCoord, locCoord, NekConstants::kGeomFactorsTol);
365 }
366 else
367 {
368 idList[i] = pFields[0]->GetExpIndex(gloCoord, locCoord,
370 }
371
372 for (int j = 0; j < locCoord.size(); ++j)
373 {
374 locCoord[j] = std::max(locCoord[j], -1.0);
375 locCoord[j] = std::min(locCoord[j], 1.0);
376 }
377
378 // If this process owns an element that contains the history point, we
379 // flag this by updating the corresponding element in procList
380 if (idList[i] != -1 && vColumnRank == 0)
381 {
382 procList[i] = vRank;
383 }
384
385 // Save a copy of the local coordinates
386 LocCoords.push_back(
387 Array<OneD, const NekDouble>(locCoord.size(), locCoord.data()));
388 }
389
390 // Clear list before populating it
391 m_historyList.clear();
392
393 // Share process ranks among all processes. If multiple processes found
394 // that they own a history point, procList will only contain the rank of
395 // the process with the largest rank
396 vComm->AllReduce(procList, LibUtilities::ReduceMax);
397 for (i = 0; i < m_historyPointsSize; ++i)
398 {
399 // If process owns a history point, it saves a reference to its global
400 // and local coordinates. Note that, if multiple processes have found
401 // that they own a history point, which can happen if the point is
402 // located on a partition boundary, only the process with the highest
403 // rank will store the point.
404 if (procList[i] == vRowRank)
405 {
406 // Note that the reference counting in the Array class prevents the
407 // data stored in GloCoords and LocCoords from being deallocated
408 // when they go out of scope, since m_historyList will increase the
409 // reference counter by 1 when it gets a reference to the data.
410 m_historyList.push_back(
411 std::make_tuple(GloCoords[i], LocCoords[i], i, idList[i]));
412 }
413 }
414
415 // Root process must have access to all coordinates for writing the file
416 if (vComm->TreatAsRankZero())
417 {
419 Array<OneD, Array<OneD, const NekDouble>>(m_historyPointsSize);
420
421 for (i = 0; i < m_historyPointsSize; ++i)
422 {
423 // Create a reference to the underlying array. Note that the data
424 // stored in GloCoords won't be deallocated because the Array class
425 // has a reference counter.
426 m_historyPoints[i] = GloCoords[i];
427 }
428 }
429
430 // Check that each history point is allocated to a process
431 if (vComm->TreatAsRankZero())
432 {
433 for (i = 0; i < m_historyPointsSize; ++i)
434 {
435 gloCoord = m_historyPoints[i];
436
437 // Write an error if no process owns history point
438 ASSERTL0(procList[i] != -1,
439 "History point " +
440 boost::lexical_cast<std::string>(gloCoord[0]) + ", " +
441 boost::lexical_cast<std::string>(gloCoord[1]) + ", " +
442 boost::lexical_cast<std::string>(gloCoord[2]) +
443 " cannot be found in the mesh.");
444 }
445
446 m_session->MatchSolverInfo("Driver", "Adaptive", m_adaptive, false);
447 }
448
450 {
451 v_Update(pFields, time);
452 }
453}
Array< OneD, Array< OneD, const NekDouble > > m_historyPoints
bool GetPoint(Array< OneD, NekDouble > gloCoord, int I)
SOLVER_UTILS_EXPORT void v_Update(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time) override
std::list< std::tuple< Array< OneD, const NekDouble >, Array< OneD, const NekDouble >, int, int > > m_historyList
std::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:55
static const NekDouble kVertexTheSameDouble
static const NekDouble kGeomFactorsTol

References ASSERTL0, GetPoint(), Nektar::NekConstants::kGeomFactorsTol, Nektar::NekConstants::kVertexTheSameDouble, m_adaptive, m_historyList, m_historyPoints, m_historyPointsSize, m_historyPointStream, m_isHomogeneous1D, m_outputPlane, m_planeIDs, Nektar::SolverUtils::Filter::m_session, Nektar::SolverUtils::Filter::m_updateOnInitialise, m_waveSpace, Nektar::LibUtilities::ReduceMax, and v_Update().

◆ v_IsTimeDependent()

bool Nektar::SolverUtils::FilterHistoryPoints::v_IsTimeDependent ( )
overrideprotectedvirtual

Implements Nektar::SolverUtils::Filter.

Definition at line 697 of file FilterHistoryPoints.cpp.

698{
699 return true;
700}

◆ v_Update()

void Nektar::SolverUtils::FilterHistoryPoints::v_Update ( const Array< OneD, const MultiRegions::ExpListSharedPtr > &  pFields,
const NekDouble time 
)
overrideprotectedvirtual

Implements Nektar::SolverUtils::Filter.

Reimplemented in Nektar::FilterCellHistoryPoints.

Definition at line 458 of file FilterHistoryPoints.cpp.

461{
462 // Only output every m_outputFrequency.
463 if ((m_index++) % m_outputFrequency)
464 {
465 return;
466 }
467
468 const int numFields = pFields.size();
469 const int coordim = pFields[0]->GetGraph()->GetSpaceDimension();
470 LibUtilities::CommSharedPtr vComm = pFields[0]->GetComm();
471 Array<OneD, NekDouble> data(m_historyPointsSize * numFields, 0.0);
472 Array<OneD, NekDouble> physvals;
473 Array<OneD, NekDouble> locCoord;
474
475 // Pull out data values field by field
476 Array<OneD, NekDouble> gloCoord(3, 0.0);
477 for (int j = 0; j < numFields; ++j)
478 {
480 {
481 Array<OneD, const unsigned int> planes = pFields[j]->GetZIDs();
482 int nPlanes = pFields[j]->GetHomogeneousBasis()->GetZ().size();
483 NekDouble lHom = pFields[j]->GetHomoLen();
484 for (auto &x : m_historyList)
485 {
486 ASSERTL0(pFields[j]->GetWaveSpace(),
487 "HistoryPoints in Homogeneous1D require that solution "
488 "is in wavespace");
489 gloCoord = std::get<0>(x);
490 locCoord = std::get<1>(x);
491 const int indx = std::get<2>(x);
492 const int expId = std::get<3>(x);
493
494 NekDouble value = 0.0;
495 NekDouble BetaT =
496 2. * M_PI * fmod(gloCoord[coordim], lHom) / lHom;
497 for (size_t n = 0; n < planes.size(); ++n)
498 {
499 if (m_waveSpace && planes[n] != m_planeIDs[indx])
500 {
501 continue;
502 }
503 physvals = pFields[j]->GetPlane(n)->UpdatePhys() +
504 pFields[j]->GetPhys_Offset(expId);
505
506 // transform elemental data if required.
507 if (pFields[j]->GetPhysState() == false)
508 {
509 pFields[j]->GetPlane(n)->GetExp(expId)->BwdTrans(
510 pFields[j]->GetPlane(n)->GetCoeffs() +
511 pFields[j]->GetCoeff_Offset(expId),
512 physvals);
513 }
514 // Interpolate data
515 NekDouble coeff =
516 pFields[j]->GetPlane(n)->GetExp(expId)->StdPhysEvaluate(
517 locCoord, physvals);
518
519 if (m_waveSpace)
520 {
521 value = coeff;
522 }
523 else
524 {
525 if (planes[n] == 0)
526 {
527 value += coeff;
528 }
529 else if (planes[n] == 1)
530 {
531 value += cos(0.5 * nPlanes * BetaT) * coeff;
532 }
533 else if (planes[n] % 2 == 0)
534 {
535 NekDouble phase = (planes[n] >> 1) * BetaT;
536 value += cos(phase) * coeff;
537 }
538 else
539 {
540 NekDouble phase = (planes[n] >> 1) * BetaT;
541 value += -sin(phase) * coeff;
542 }
543 }
544 }
545 // store data
546 data[indx * numFields + j] = value;
547 }
548 }
549 else
550 {
551 for (auto &x : m_historyList)
552 {
553 locCoord = std::get<1>(x);
554 const int indx = std::get<2>(x);
555 const int expId = std::get<3>(x);
556
557 physvals = pFields[j]->UpdatePhys() +
558 pFields[j]->GetPhys_Offset(expId);
559
560 // transform elemental data if required.
561 if (pFields[j]->GetPhysState() == false)
562 {
563 pFields[j]->GetExp(expId)->BwdTrans(
564 pFields[j]->GetCoeffs() +
565 pFields[j]->GetCoeff_Offset(expId),
566 physvals);
567 }
568
569 // interpolate point
570 data[indx * numFields + j] =
571 pFields[j]->GetExp(expId)->StdPhysEvaluate(locCoord,
572 physvals);
573 }
574 }
575 }
576
577 // Exchange history data
578 // This could be improved to reduce communication but works for now
579 vComm->AllReduce(data, LibUtilities::ReduceSum);
580
581 // TODO: Why not only call this routine if we are rank 0?
582 v_WriteData(vComm->GetRank(), data, numFields, time);
583}
SOLVER_UTILS_EXPORT void v_WriteData(const int &rank, const Array< OneD, NekDouble > &data, const int &numFields, const NekDouble &time)

References ASSERTL0, m_historyList, m_historyPointsSize, m_index, m_isHomogeneous1D, m_outputFrequency, m_planeIDs, m_waveSpace, Nektar::LibUtilities::ReduceSum, and v_WriteData().

Referenced by v_Initialise().

◆ v_WriteData()

void Nektar::SolverUtils::FilterHistoryPoints::v_WriteData ( const int &  rank,
const Array< OneD, NekDouble > &  data,
const int &  numFields,
const NekDouble time 
)
protected

Definition at line 585 of file FilterHistoryPoints.cpp.

589{
590 // Only the root process writes out history data
591 if (rank == 0)
592 {
593 Array<OneD, NekDouble> gloCoord(3, 0.0);
594 if (!m_outputOneFile || m_index == 1)
595 {
596 std::stringstream vTmpFilename;
597 std::string vOutputFilename;
598 // get the file extension
599 std::string ext = fs::path(m_outputFile).extension().string();
600 if (m_outputOneFile)
601 {
602 vTmpFilename << m_outputFile;
603 }
604 else
605 {
606 vTmpFilename
607 << fs::path(m_outputFile).replace_extension("").string()
608 << "_" << m_outputIndex << ext;
609 }
610 // back up the file if already exists and backup switch is turned on
611 vOutputFilename = Filter::SetupOutput(ext, vTmpFilename.str());
612
614 if (m_adaptive)
615 {
616 m_outputStream.open(vOutputFilename.c_str(), ofstream::app);
617 }
618 else
619 {
620 m_outputStream.open(vOutputFilename.c_str());
621 }
622 m_outputStream << "# History data for variables (:";
623
624 for (int i = 0; i < numFields; ++i)
625 {
626 m_outputStream << m_session->GetVariable(i) << ",";
627 }
628
630 {
631 m_outputStream << ") at points:" << endl;
632 }
633 else
634 {
635 m_outputStream << ") at points:" << endl;
636 }
637
638 for (int i = 0; i < m_historyPoints.size(); ++i)
639 {
640 gloCoord = m_historyPoints[i];
641
642 m_outputStream << "# " << boost::format("%6.0f") % i;
643 m_outputStream << " " << boost::format("%25.19e") % gloCoord[0];
644 m_outputStream << " " << boost::format("%25.19e") % gloCoord[1];
645 m_outputStream << " " << boost::format("%25.19e") % gloCoord[2];
646 m_outputStream << endl;
647 }
648
650 {
651 if (m_waveSpace)
652 {
653 m_outputStream << "# (in Wavespace)" << endl;
654 }
655 }
656 }
657
658 // Write data values point by point
659 for (int k = 0; k < m_historyPoints.size(); ++k)
660 {
661 m_outputStream << boost::format("%25.19e") % time;
662 for (int j = 0; j < numFields; ++j)
663 {
665 << " "
666 << boost::format("%25.19e") % data[k * numFields + j];
667 }
668 m_outputStream << endl;
669 }
670
671 if (!m_outputOneFile)
672 {
673 m_outputStream.close();
674 }
675 }
676}

References CellMLToNektar.pycml::format, m_adaptive, m_historyPoints, m_index, m_isHomogeneous1D, m_outputFile, m_outputIndex, m_outputOneFile, m_outputStream, Nektar::SolverUtils::Filter::m_session, m_waveSpace, and Nektar::SolverUtils::Filter::SetupOutput().

Referenced by v_Update().

Friends And Related Function Documentation

◆ MemoryManager< FilterHistoryPoints >

friend class MemoryManager< FilterHistoryPoints >
friend

Definition at line 1 of file FilterHistoryPoints.h.

Member Data Documentation

◆ className

std::string Nektar::SolverUtils::FilterHistoryPoints::className
static
Initial value:
=
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
static FilterSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const std::shared_ptr< EquationSystem > &pEquation, const std::map< std::string, std::string > &pParams)
Creates an instance of this class.
FilterFactory & GetFilterFactory()

Name of the class.

Definition at line 63 of file FilterHistoryPoints.h.

◆ m_adaptive

bool Nektar::SolverUtils::FilterHistoryPoints::m_adaptive
protected

◆ m_historyList

std::list<std::tuple<Array<OneD, const NekDouble>, Array<OneD, const NekDouble>, int, int> > Nektar::SolverUtils::FilterHistoryPoints::m_historyList
protected

◆ m_historyPoints

Array<OneD, Array<OneD, const NekDouble> > Nektar::SolverUtils::FilterHistoryPoints::m_historyPoints
protected
Initial value:
=
Array<OneD, Array<OneD, const NekDouble>>(0)

Definition at line 88 of file FilterHistoryPoints.h.

Referenced by v_Initialise(), Nektar::FilterCellHistoryPoints::v_Update(), v_WriteData(), and Nektar::FilterCellHistoryPoints::v_WriteData().

◆ m_historyPointsSize

size_t Nektar::SolverUtils::FilterHistoryPoints::m_historyPointsSize = 0
protected

Definition at line 90 of file FilterHistoryPoints.h.

Referenced by v_Initialise(), and v_Update().

◆ m_historyPointStream

std::stringstream Nektar::SolverUtils::FilterHistoryPoints::m_historyPointStream
protected

Definition at line 100 of file FilterHistoryPoints.h.

Referenced by FilterHistoryPoints(), GetPoint(), and v_Initialise().

◆ m_index

unsigned int Nektar::SolverUtils::FilterHistoryPoints::m_index = 0
protected

◆ m_isHomogeneous1D

bool Nektar::SolverUtils::FilterHistoryPoints::m_isHomogeneous1D
protected

◆ m_outputFile

std::string Nektar::SolverUtils::FilterHistoryPoints::m_outputFile
protected

◆ m_outputFrequency

unsigned int Nektar::SolverUtils::FilterHistoryPoints::m_outputFrequency
protected

◆ m_outputIndex

unsigned int Nektar::SolverUtils::FilterHistoryPoints::m_outputIndex = 0
protected

◆ m_outputOneFile

bool Nektar::SolverUtils::FilterHistoryPoints::m_outputOneFile
protected

◆ m_outputPlane

int Nektar::SolverUtils::FilterHistoryPoints::m_outputPlane
protected

plane to take history point from if using a homogeneous1D expansion

Definition at line 94 of file FilterHistoryPoints.h.

Referenced by FilterHistoryPoints(), v_Initialise(), and Nektar::FilterCellHistoryPoints::v_Update().

◆ m_outputStream

std::ofstream Nektar::SolverUtils::FilterHistoryPoints::m_outputStream
protected

◆ m_planeIDs

std::vector<int> Nektar::SolverUtils::FilterHistoryPoints::m_planeIDs
protected

Definition at line 95 of file FilterHistoryPoints.h.

Referenced by v_Initialise(), and v_Update().

◆ m_pointDatMap

std::map<LibUtilities::PtsType, Array<OneD, NekDouble> > Nektar::SolverUtils::FilterHistoryPoints::m_pointDatMap
protected

Definition at line 110 of file FilterHistoryPoints.h.

Referenced by FilterHistoryPoints(), and GetPoint().

◆ m_pointNumMap

std::map<LibUtilities::PtsType, Array<OneD, int> > Nektar::SolverUtils::FilterHistoryPoints::m_pointNumMap
protected

Definition at line 111 of file FilterHistoryPoints.h.

Referenced by FilterHistoryPoints(), and GetPoint().

◆ m_waveSpace

bool Nektar::SolverUtils::FilterHistoryPoints::m_waveSpace
protected