diff --git a/CHANGELOG.md b/CHANGELOG.md
index ede9136c5c91bb0fac3230a4733d54f2ad1cbf7c..b995c7bbd764743583f97c333798115da06d1d6d 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -39,7 +39,9 @@ v5.8.0
- Updated the User-guide with additional inofrmation for outflow BC, addressing the issue #103 (!1988)
- Updated the User-guide with additional inofrmation for outflow BC, addressing the issue #103 (!1990)
-
+**PulseWaveSolver**
+- Added ability to output history points (and other filters) (!2000)
+
v5.7.0
-----
**Library**
diff --git a/library/LibUtilities/BasicUtils/SessionReader.cpp b/library/LibUtilities/BasicUtils/SessionReader.cpp
index 882d31fe922c8a13992607cc10e1f3ac221c56e1..c40a37ee9477f020c38aa737a056ae4d2b6d6b5b 100644
--- a/library/LibUtilities/BasicUtils/SessionReader.cpp
+++ b/library/LibUtilities/BasicUtils/SessionReader.cpp
@@ -2409,6 +2409,15 @@ void SessionReader::ReadFilters(TiXmlElement *filters)
"Missing attribute 'TYPE' for filter.");
std::string typeStr = filter->Attribute("TYPE");
+ int domainID = -1;
+ if (filter->Attribute("DOMAIN"))
+ {
+ // domainID = int(filter->Attribute("DOMAIN"));
+ int err = filter->QueryIntAttribute("DOMAIN", &domainID);
+ ASSERTL0(err == TIXML_SUCCESS,
+ "Unable to read attribute DOMAIN in filter.");
+ }
+
std::map<std::string, std::string> vParams;
TiXmlElement *element = filter;
@@ -2429,8 +2438,11 @@ void SessionReader::ReadFilters(TiXmlElement *filters)
param = param->NextSiblingElement("PARAM");
}
- m_filters.push_back(
- std::pair<std::string, FilterParams>(typeStr, vParams));
+ FilterDefinition filterDef;
+ filterDef.domain = domainID;
+ filterDef.name = typeStr;
+ filterDef.params = vParams;
+ m_filters.push_back(filterDef);
filter = filter->NextSiblingElement("FILTER");
}
diff --git a/library/LibUtilities/BasicUtils/SessionReader.h b/library/LibUtilities/BasicUtils/SessionReader.h
index 86071be273bc70fc48c350366d2cd48f6c487a98..e016103e4126451daf9c7f15a02d248e7a42a8ad 100644
--- a/library/LibUtilities/BasicUtils/SessionReader.h
+++ b/library/LibUtilities/BasicUtils/SessionReader.h
@@ -59,7 +59,13 @@ typedef std::map<std::string, std::string> GeometricInfoMap;
typedef std::vector<std::string> VariableList;
typedef std::map<std::string, std::string> TagMap;
typedef std::map<std::string, std::string> FilterParams;
-typedef std::vector<std::pair<std::string, FilterParams>> FilterMap;
+struct FilterDefinition
+{
+ std::string name;
+ int domain;
+ FilterParams params;
+};
+typedef std::vector<FilterDefinition> FilterMap;
struct CmdLineArg
{
diff --git a/library/SolverUtils/Filters/Filter.h b/library/SolverUtils/Filters/Filter.h
index a25eefac8e21e75031a5281b842133e20d1fd857..c6446a3fdfb7d2ea7848585d592faffa6952a142 100644
--- a/library/SolverUtils/Filters/Filter.h
+++ b/library/SolverUtils/Filters/Filter.h
@@ -84,9 +84,15 @@ public:
SOLVER_UTILS_EXPORT inline std::string SetupOutput(
const std::string ext, const std::string inname);
+ SOLVER_UTILS_EXPORT inline void SetUpdateOnInitialise(bool flag)
+ {
+ m_updateOnInitialise = flag;
+ }
+
protected:
LibUtilities::SessionReaderSharedPtr m_session;
const std::weak_ptr<EquationSystem> m_equ;
+ bool m_updateOnInitialise = true;
virtual void v_Initialise(
const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields,
diff --git a/library/SolverUtils/Filters/FilterAeroForces.cpp b/library/SolverUtils/Filters/FilterAeroForces.cpp
index e3a2cd1fcd53da87a872988cdee960587db05e9a..57ac9eee6dd1f9b4c002e849c780f771dcdbe555 100644
--- a/library/SolverUtils/Filters/FilterAeroForces.cpp
+++ b/library/SolverUtils/Filters/FilterAeroForces.cpp
@@ -220,7 +220,7 @@ FilterAeroForces::~FilterAeroForces()
*/
void FilterAeroForces::v_Initialise(
const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields,
- const NekDouble &time)
+ [[maybe_unused]] const NekDouble &time)
{
// Load mapping
m_mapping = GlobalMapping::Mapping::Load(m_session, pFields);
@@ -418,7 +418,11 @@ void FilterAeroForces::v_Initialise(
m_lastTime = -1;
m_index = 0;
- v_Update(pFields, time);
+
+ if (m_updateOnInitialise)
+ {
+ v_Update(pFields, time);
+ }
}
/**
diff --git a/library/SolverUtils/Filters/FilterCheckpoint.cpp b/library/SolverUtils/Filters/FilterCheckpoint.cpp
index 3f30aa710b2be1c67942ccfb23800709dd2327c4..e6c9073b7ddc3ddca0363d2f6daa6924ee6ab69e 100644
--- a/library/SolverUtils/Filters/FilterCheckpoint.cpp
+++ b/library/SolverUtils/Filters/FilterCheckpoint.cpp
@@ -75,12 +75,17 @@ FilterCheckpoint::~FilterCheckpoint()
}
void FilterCheckpoint::v_Initialise(
- const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields,
- const NekDouble &time)
+ [[maybe_unused]] const Array<OneD, const MultiRegions::ExpListSharedPtr>
+ &pFields,
+ [[maybe_unused]] const NekDouble &time)
{
m_index = 0;
m_outputIndex = 0;
- v_Update(pFields, time);
+
+ if (m_updateOnInitialise)
+ {
+ v_Update(pFields, time);
+ }
}
void FilterCheckpoint::v_Update(
diff --git a/library/SolverUtils/Filters/FilterEnergy.cpp b/library/SolverUtils/Filters/FilterEnergy.cpp
index 3188ac4c27b97fae77980adeea10c1b428dc5c0d..0979b75e70da50c10268f705189606e8a813254a 100644
--- a/library/SolverUtils/Filters/FilterEnergy.cpp
+++ b/library/SolverUtils/Filters/FilterEnergy.cpp
@@ -80,7 +80,7 @@ FilterEnergy::~FilterEnergy()
void FilterEnergy::v_Initialise(
const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields,
- const NekDouble &time)
+ [[maybe_unused]] const NekDouble &time)
{
m_index = -1;
MultiRegions::ExpListSharedPtr areaField;
@@ -117,9 +117,10 @@ void FilterEnergy::v_Initialise(
{
m_area *= m_homogeneousLength;
}
-
- // Output values at initial time.
- v_Update(pFields, time);
+ if (m_updateOnInitialise)
+ {
+ v_Update(pFields, time);
+ }
}
void FilterEnergy::v_Update(
diff --git a/library/SolverUtils/Filters/FilterFieldConvert.cpp b/library/SolverUtils/Filters/FilterFieldConvert.cpp
index 9aea7af4832823ce5fbc8a9567604842182b1550..b503e5c43d8e8693edb907dd2c1c611cf0e046c6 100644
--- a/library/SolverUtils/Filters/FilterFieldConvert.cpp
+++ b/library/SolverUtils/Filters/FilterFieldConvert.cpp
@@ -439,6 +439,7 @@ void FilterFieldConvert::v_Update(
const NekDouble &time)
{
m_index++;
+
if (m_index % m_sampleFrequency > 0 ||
(time - m_outputStartTime) < -1.0e-07)
{
diff --git a/library/SolverUtils/Filters/FilterHistoryPoints.cpp b/library/SolverUtils/Filters/FilterHistoryPoints.cpp
index ff199b3c018bef7be3f922a7763e3918d10c53b9..6dfbeb80a8610c809c46aa36cc7dd6e3dc8a5a87 100644
--- a/library/SolverUtils/Filters/FilterHistoryPoints.cpp
+++ b/library/SolverUtils/Filters/FilterHistoryPoints.cpp
@@ -288,7 +288,7 @@ bool FilterHistoryPoints::GetPoint(Array<OneD, NekDouble> gloCoord, int I)
*/
void FilterHistoryPoints::v_Initialise(
const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields,
- const NekDouble &time)
+ [[maybe_unused]] const NekDouble &time)
{
ASSERTL0(!m_historyPointStream.fail(), "No history points in stream.");
@@ -445,7 +445,11 @@ void FilterHistoryPoints::v_Initialise(
m_session->MatchSolverInfo("Driver", "Adaptive", m_adaptive, false);
}
- v_Update(pFields, time);
+
+ if (m_updateOnInitialise)
+ {
+ v_Update(pFields, time);
+ }
}
/**
@@ -667,7 +671,6 @@ void FilterHistoryPoints::v_WriteData(const int &rank,
if (!m_outputOneFile)
{
m_outputStream.close();
- ;
}
}
}
@@ -681,7 +684,10 @@ void FilterHistoryPoints::v_Finalise(
{
if (pFields[0]->GetComm()->GetRank() == 0 && m_outputOneFile)
{
- m_outputStream.close();
+ if (m_outputStream.is_open())
+ {
+ m_outputStream.close();
+ }
}
}
diff --git a/library/SolverUtils/Filters/FilterIntegral.cpp b/library/SolverUtils/Filters/FilterIntegral.cpp
index 83b493c7e846bee8bf0ee7d37d2dbf0c3d0a3e7f..bd52d3c33af2936d0a59da9f1d2de64ffc99136f 100644
--- a/library/SolverUtils/Filters/FilterIntegral.cpp
+++ b/library/SolverUtils/Filters/FilterIntegral.cpp
@@ -143,7 +143,7 @@ FilterIntegral::FilterIntegral(
*/
void FilterIntegral::v_Initialise(
const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields,
- const NekDouble &time)
+ [[maybe_unused]] const NekDouble &time)
{
// Create map from element ID -> expansion list ID
@@ -266,8 +266,10 @@ void FilterIntegral::v_Initialise(
m_compExpMap[i] = tmpCompExp;
}
-
- v_Update(pFields, time);
+ if (m_updateOnInitialise)
+ {
+ v_Update(pFields, time);
+ }
}
/**
diff --git a/library/SolverUtils/Filters/FilterLagrangianPoints.cpp b/library/SolverUtils/Filters/FilterLagrangianPoints.cpp
index 49ac990d520d5829e27d763c985dcbf01f02fdc6..449c7aaa177215142145dca443cd09d9f48d751c 100644
--- a/library/SolverUtils/Filters/FilterLagrangianPoints.cpp
+++ b/library/SolverUtils/Filters/FilterLagrangianPoints.cpp
@@ -659,7 +659,10 @@ void FilterLagrangianPoints::v_Initialise(
m_ofstreamSamplePoints << endl;
}
SetUpCommInfo();
- v_Update(pFields, time);
+ if (m_updateOnInitialise)
+ {
+ v_Update(pFields, time);
+ }
}
void FilterLagrangianPoints::v_ModifyVelocity(
diff --git a/library/SolverUtils/Filters/FilterMean.cpp b/library/SolverUtils/Filters/FilterMean.cpp
index 02cd9bca6b5eedec96ff252d797f9d3192c65461..5d401d6595720710b7215c79e847f9bebd3b58c3 100644
--- a/library/SolverUtils/Filters/FilterMean.cpp
+++ b/library/SolverUtils/Filters/FilterMean.cpp
@@ -69,7 +69,7 @@ FilterMean::~FilterMean()
void FilterMean::v_Initialise(
const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields,
- const NekDouble &time)
+ [[maybe_unused]] const NekDouble &time)
{
MultiRegions::ExpListSharedPtr areaField;
int spacedim = 2;
@@ -132,7 +132,10 @@ void FilterMean::v_Initialise(
// Output values at initial time.
m_index = 0;
- v_Update(pFields, time);
+ if (m_updateOnInitialise)
+ {
+ v_Update(pFields, time);
+ }
}
void FilterMean::v_Update(
diff --git a/library/SolverUtils/Filters/FilterModalEnergy.cpp b/library/SolverUtils/Filters/FilterModalEnergy.cpp
index fb2e08fbd1721bae230235b77bdda19229171cfc..dc142744ff5a17ce81f44e126eaed8d7d75f0dd5 100644
--- a/library/SolverUtils/Filters/FilterModalEnergy.cpp
+++ b/library/SolverUtils/Filters/FilterModalEnergy.cpp
@@ -109,7 +109,7 @@ FilterModalEnergy::~FilterModalEnergy()
*/
void FilterModalEnergy::v_Initialise(
const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields,
- const NekDouble &time)
+ [[maybe_unused]] const NekDouble &time)
{
LibUtilities::CommSharedPtr vComm = pFields[0]->GetComm();
@@ -139,7 +139,10 @@ void FilterModalEnergy::v_Initialise(
}
m_index = 0;
- v_Update(pFields, time);
+ if (m_updateOnInitialise)
+ {
+ v_Update(pFields, time);
+ }
}
/**
diff --git a/library/SolverUtils/UnsteadySystem.cpp b/library/SolverUtils/UnsteadySystem.cpp
index f0bde5e9abe6d33356047cf213d3ad09ef237755..b48f541e29c90c2dfdd530b4435612dc8cf20d20 100644
--- a/library/SolverUtils/UnsteadySystem.cpp
+++ b/library/SolverUtils/UnsteadySystem.cpp
@@ -165,8 +165,8 @@ void UnsteadySystem::v_InitObject(bool DeclareField)
for (auto &x : m_session->GetFilters())
{
m_filters.push_back(make_pair(
- x.first, GetFilterFactory().CreateInstance(
- x.first, m_session, shared_from_this(), x.second)));
+ x.name, GetFilterFactory().CreateInstance(
+ x.name, m_session, shared_from_this(), x.params)));
}
}
diff --git a/solvers/IncNavierStokesSolver/EquationSystems/SmoothedProfileMethod.cpp b/solvers/IncNavierStokesSolver/EquationSystems/SmoothedProfileMethod.cpp
index ac34856ad133f52ddb8595a5be9b24057300b6df..0852df57aee5f6a91d1037de0734d5d06adbcfe1 100644
--- a/solvers/IncNavierStokesSolver/EquationSystems/SmoothedProfileMethod.cpp
+++ b/solvers/IncNavierStokesSolver/EquationSystems/SmoothedProfileMethod.cpp
@@ -217,7 +217,7 @@ void SmoothedProfileMethod::v_InitObject(bool DeclareField)
m_forcesFilter = -1;
for (size_t i = 0; i < m_session->GetFilters().size(); ++i)
{
- if (boost::iequals(m_session->GetFilters()[i].first, "AeroForcesSPM"))
+ if (boost::iequals(m_session->GetFilters()[i].name, "AeroForcesSPM"))
{
m_forcesFilter = i;
break;
diff --git a/solvers/PulseWaveSolver/EquationSystems/PulseWavePropagation.cpp b/solvers/PulseWaveSolver/EquationSystems/PulseWavePropagation.cpp
index 50c82acd2a457b90ff0f88ec80cacc86c40f703d..1baacc826e7f9fabf0a08b2884ed646a8cf5b6e0 100644
--- a/solvers/PulseWaveSolver/EquationSystems/PulseWavePropagation.cpp
+++ b/solvers/PulseWaveSolver/EquationSystems/PulseWavePropagation.cpp
@@ -177,7 +177,7 @@ void PulseWavePropagation::DoOdeRhs(
timer.Start();
for (i = 0; i < m_nVariables; ++i)
{
- physarray[i] = inarray[i] + cnt;
+ physarray[i] = inarray[i] + cnt; // note this is doing a hidden copy
out[i] = outarray[i] + cnt;
}
@@ -319,13 +319,13 @@ void PulseWavePropagation::GetFluxVector(
LibUtilities::Timer timer;
+ timer.Start();
for (size_t j = 0; j < nq; ++j)
{
- timer.Start();
flux[0][0][j] = physfield[0][j] * physfield[1][j];
- timer.Stop();
- timer.AccumulateRegion("PulseWavePropagation:GetFluxVector-flux", 3);
}
+ timer.Stop();
+ timer.AccumulateRegion("PulseWavePropagation:GetFluxVector-flux", 3);
// d/dx of AU, for the viscoelastic tube law and extra fields
m_fields[0]->PhysDeriv(flux[0][0], dAUdx);
diff --git a/solvers/PulseWaveSolver/EquationSystems/PulseWaveSystem.cpp b/solvers/PulseWaveSolver/EquationSystems/PulseWaveSystem.cpp
index d4e5840e6159617043cd622193a248d19155818d..289edde603e600e59577ac793fa065a62026d5f9 100644
--- a/solvers/PulseWaveSolver/EquationSystems/PulseWaveSystem.cpp
+++ b/solvers/PulseWaveSolver/EquationSystems/PulseWaveSystem.cpp
@@ -132,28 +132,62 @@ void PulseWaveSystem::v_InitObject(bool DeclareField)
SpatialDomains::BoundaryConditions Allbcs(m_session, m_graph);
- // Set up domains and put geometry to be only one space dimension.
- size_t cnt = 0;
- bool SetToOneSpaceDimension = true;
-
- if (m_session->DefinesCmdLineArgument("SetToOneSpaceDimension"))
- {
- std::string cmdline = m_session->GetCmdLineArgument<std::string>(
- "SetToOneSpaceDimension");
- if (boost::to_upper_copy(cmdline) == "FALSE")
- {
- SetToOneSpaceDimension = false;
- }
- }
-
// get parallel communicators as required
map<int, LibUtilities::CommSharedPtr> domComm;
GetCommArray(domComm);
SetUpDomainInterfaceBCs(Allbcs, domComm);
+ // set up a list of the domain to filter ids
+ auto filterMap = m_session->GetFilters();
+ unsigned int fcnt = 0;
+ for (auto &x : filterMap)
+ {
+ ASSERTL0(x.domain != -1, "Filter needs DOMAIN attribute when "
+ "applied to PulseWaveSolver.")
+ m_domainToFilterIDs[x.domain].push_back(fcnt);
+ ++fcnt;
+ }
+
+ // Set up domains and put geometry to be only one space dimension.
+ size_t cnt = 0;
+ // currently I believe we have to set to one space dimension on
+ // constuction since the advection terms are set to differentiate
+ // with respect to x (or the lcoal ordinate) not as an arc length
+ // derivative. So will not get consistent answers wtihout this
+ bool SetToOneSpaceDimension = true;
+
for (auto &d : m_domOrder)
{
+ // need to check if we have filters on this domain and if so
+ // set them up without SetToOneSpaceDimension if so we can
+ // initialisatio them (i.e. history points)
+ if (m_domainToFilterIDs.count(d))
+ {
+ // make temporary expansion list to not settig to one space
+ // dimension
+ Array<OneD, MultiRegions::ExpListSharedPtr> filterFields(
+ m_nVariables);
+ for (size_t j = 0; j < m_nVariables; ++j)
+ {
+ filterFields[j] = MemoryManager<MultiRegions::DisContField>::
+ AllocateSharedPtr(m_session, m_graph, m_domain[d], Allbcs,
+ m_session->GetVariable(j), domComm[d],
+ false);
+ }
+
+ // intiialise filters with temporary expansion lsits
+ for (auto &x : m_domainToFilterIDs[d])
+ {
+ m_filters[x].second->SetUpdateOnInitialise(false);
+ m_filters[x].second->Initialise(filterFields, m_time);
+ m_filters[x].second->SetUpdateOnInitialise(true);
+
+ // save vessel start id for update
+ m_filterToVesselID[x] = cnt;
+ }
+ }
+
for (size_t j = 0; j < m_nVariables; ++j)
{
m_vessels[cnt++] =
@@ -295,26 +329,23 @@ void PulseWaveSystem::v_InitObject(bool DeclareField)
m_alpha_trace[omega] = Array<OneD, NekDouble>(nqTrace);
m_fields[0]->ExtractTracePhys(m_alpha[omega], m_alpha_trace[omega]);
- if (SetToOneSpaceDimension)
- {
- m_trace_fwd_normal[omega] = Array<OneD, NekDouble>(nqTrace, 0.0);
+ m_trace_fwd_normal[omega] = Array<OneD, NekDouble>(nqTrace, 0.0);
- MultiRegions::ExpListSharedPtr trace = m_fields[0]->GetTrace();
- int nelmt_trace = trace->GetExpSize();
+ MultiRegions::ExpListSharedPtr trace = m_fields[0]->GetTrace();
+ int nelmt_trace = trace->GetExpSize();
- Array<OneD, Array<OneD, NekDouble>> normals(nelmt_trace);
+ Array<OneD, Array<OneD, NekDouble>> normals(nelmt_trace);
- for (int i = 0; i < nelmt_trace; ++i)
- {
- normals[i] = m_trace_fwd_normal[omega] + i;
- }
+ for (int i = 0; i < nelmt_trace; ++i)
+ {
+ normals[i] = m_trace_fwd_normal[omega] + i;
+ }
- // need to set to 1 for consistency since boundary
- // conditions may not have coordim = 1
- trace->GetExp(0)->GetGeom()->SetCoordim(1);
+ // need to set to 1 for consistency since boundary
+ // conditions may not have coordim = 1
+ trace->GetExp(0)->GetGeom()->SetCoordim(1);
- trace->GetNormals(normals);
- }
+ trace->GetNormals(normals);
omega++;
}
@@ -458,8 +489,8 @@ void PulseWaveSystem::GetCommArray(
while (search)
{
size_t dorank = 0; // search index over processors
- set<int>
- proclist; // has proc been identified for a comm at this level
+ set<int> proclist; // has proc been identified for a comm at
+ // this level
int flag = 100; // colour for communicators in this search level
while (dorank < nprocs)
@@ -683,7 +714,7 @@ void PulseWaveSystem::GetCommArray(
Array<OneD, NekDouble> shareddom(totShared, -1.0);
cnt = 0; // collect a list of domains that each shared commm is
- // attached to
+ // attached to
for (auto &s : SharedProc)
{
shareddom[numShared[rank] + cnt] = s.first;
@@ -736,7 +767,8 @@ void PulseWaveSystem::GetCommArray(
m_comm->AllReduce(maxoffset, LibUtilities::ReduceMax);
maxoffset++;
- // make a map relating the order of the SharedProc to the domain id
+ // make a map relating the order of the SharedProc to the domain
+ // id
map<int, int> ShrToDom;
cnt = 0;
for (auto &s : SharedProc)
@@ -1176,8 +1208,8 @@ void PulseWaveSystem::SetUpDomainInterfaceBCs(
MemoryManager<
SpatialDomains::BoundaryRegion>::AllocateSharedPtr());
- // Set up Composite (GemetryVector) to contain vertex and put into
- // bRegion
+ // Set up Composite (GemetryVector) to contain vertex and put
+ // into bRegion
SpatialDomains::CompositeSharedPtr gvec =
MemoryManager<SpatialDomains::Composite>::AllocateSharedPtr();
gvec->m_geomVec.push_back(b.second);
@@ -1226,7 +1258,7 @@ void PulseWaveSystem::v_DoInitialise(
int omega = 0;
for (auto &d : m_domOrder)
{
- for (size_t i = 0; i < 2; ++i)
+ for (size_t i = 0; i < m_nVariables; ++i)
{
m_fields[i] = m_vessels[m_nVariables * omega + i];
}
@@ -1240,11 +1272,29 @@ void PulseWaveSystem::v_DoInitialise(
omega++;
}
- // Reset 2 variables to first vessels
- for (size_t i = 0; i < 2; ++i)
+ // Reset variables to first vessels
+ for (size_t i = 0; i < m_nVariables; ++i)
{
m_fields[i] = m_vessels[i];
}
+
+ // Update filters.
+ for (auto &d : m_domainToFilterIDs)
+ {
+ Array<OneD, MultiRegions::ExpListSharedPtr> filterFields(m_nVariables);
+
+ int vesselID = m_filterToVesselID[d.second[0]];
+ for (size_t j = 0; j < m_nVariables; ++j)
+ {
+ filterFields[j] = m_vessels[vesselID + j];
+ }
+
+ // loop over filters in domain
+ for (auto &f : d.second)
+ {
+ m_filters[f].second->Update(filterFields, m_time);
+ }
+ }
}
/**
@@ -1266,7 +1316,6 @@ void PulseWaveSystem::v_DoInitialise(
*/
void PulseWaveSystem::v_DoSolve()
{
- // NekDouble IntegrationTime = 0.0;
size_t i;
int n;
int nchk = 1;
@@ -1291,8 +1340,6 @@ void PulseWaveSystem::v_DoSolve()
time_v_IntStep.Stop();
time_v_IntStep.AccumulateRegion("PulseWaveSystem::TimeIntegrationStep",
0);
- // IntegrationTime += timer.TimePerTest(1);
-
// Write out status information.
if (m_infosteps && !((n + 1) % m_infosteps) &&
m_session->GetComm()->GetRank() == 0)
@@ -1301,38 +1348,72 @@ void PulseWaveSystem::v_DoSolve()
<< "\t Time-step: " << m_timestep << "\t" << endl;
}
+ // Copy Array To Vessel Phys Fields
+ for (size_t i = 0; i < m_nVariables; ++i)
+ {
+ Vmath::Vcopy(fields[i].size(), fields[i], 1,
+ m_vessels[i]->UpdatePhys(), 1);
+ }
+
// Transform data if needed
if (m_checksteps && !((n + 1) % m_checksteps))
{
for (i = 0; i < m_nVariables; ++i)
{
- size_t cnt = 0;
for (size_t omega = 0; omega < m_nDomains; omega++)
{
- m_vessels[omega * m_nVariables + i]->FwdTrans(
- fields[i] + cnt,
- m_vessels[omega * m_nVariables + i]->UpdateCoeffs());
- cnt += m_vessels[omega * m_nVariables + i]->GetTotPoints();
+ unsigned int offset = omega * m_nVariables + i;
+ m_vessels[offset]->FwdTrans(
+ fields[i] + m_fieldPhysOffset[omega],
+ m_vessels[offset]->UpdateCoeffs());
}
}
CheckPoint_Output(nchk++);
}
+ // Update filters.
+ for (auto &d : m_domainToFilterIDs)
+ {
+ Array<OneD, MultiRegions::ExpListSharedPtr> filterFields(
+ m_nVariables);
+ int vesselID = m_filterToVesselID[d.second[0]];
+ for (size_t j = 0; j < m_nVariables; ++j)
+ {
+ filterFields[j] = m_vessels[vesselID + j];
+ }
+
+ // loop over filters in domain
+ for (auto &f : d.second)
+ {
+ m_filters[f].second->Update(filterFields, m_time);
+ }
+ }
+
} // end of timeintegration
- // Copy Array To Vessel Phys Fields
- for (size_t i = 0; i < m_nVariables; ++i)
+ /** if(m_session->GetComm()->GetRank() == 0)
+ {
+ cout << "Time-integration timing: " << IntegrationTime << " s" <<
+ endl
+ << endl;
+ } **/
+
+ // Finalise filters.
+ for (auto &d : m_domainToFilterIDs)
{
- Vmath::Vcopy(fields[i].size(), fields[i], 1, m_vessels[i]->UpdatePhys(),
- 1);
- }
+ Array<OneD, MultiRegions::ExpListSharedPtr> filterFields(m_nVariables);
+ int vesselID = m_filterToVesselID[d.second[0]];
+ for (size_t j = 0; j < m_nVariables; ++j)
+ {
+ filterFields[j] = m_vessels[vesselID + j];
+ }
- /** if(m_session->GetComm()->GetRank() == 0)
- {
- cout << "Time-integration timing: " << IntegrationTime << " s" <<
- endl
- << endl;
- } **/
+ // loop over filters in domain
+ for (auto &f : d.second)
+ {
+ m_filters[f].second->Finalise(filterFields, m_time);
+ }
+ }
}
void PulseWaveSystem::FillDataFromInterfacePoint(
@@ -1552,15 +1633,18 @@ void PulseWaveSystem::BifurcationRiemann(Array<OneD, NekDouble> &Au,
/*
* We solve the six constraint equations via a multivariate Newton
* iteration. Equations are:
- * 1. Forward characteristic: W1(A_L, U_L) = W1(Au_L, Uu_L)
- * 2. Backward characteristic 1: W2(A_R1, U_R1) = W2(Au_R1, Uu_R1)
- * 3. Backward characteristic 2: W2(A_R2, U_R2) = W2(Au_R2, Uu_R2)
- * 4. Conservation of mass: Au_L * Uu_L = Au_R1 * Uu_R1 +
- * Au_R2 * Uu_R2
- * 5. Continuity of total pressure 1: rho * Uu_L * Uu_L / 2 + p(Au_L)
- * = rho * Uu_R1 * Uu_R1 / 2 + p(Au_R1)
- * 6. Continuity of total pressure 2: rho * Uu_L * Uu_L / 2 + p(Au_L)
- * = rho * Uu_R2 * Uu_R2 / 2 + p(Au_R2)
+ * 1. Forward characteristic: W1(A_L, U_L) = W1(Au_L,
+ * Uu_L)
+ * 2. Backward characteristic 1: W2(A_R1, U_R1) = W2(Au_R1,
+ * Uu_R1)
+ * 3. Backward characteristic 2: W2(A_R2, U_R2) = W2(Au_R2,
+ * Uu_R2)
+ * 4. Conservation of mass: Au_L * Uu_L = Au_R1 *
+ * Uu_R1 + Au_R2 * Uu_R2
+ * 5. Continuity of total pressure 1: rho * Uu_L * Uu_L / 2 +
+ * p(Au_L) = rho * Uu_R1 * Uu_R1 / 2 + p(Au_R1)
+ * 6. Continuity of total pressure 2: rho * Uu_L * Uu_L / 2 +
+ * p(Au_L) = rho * Uu_R2 * Uu_R2 / 2 + p(Au_R2)
*/
for (size_t i = 0; i < 3; ++i)
{
@@ -1662,17 +1746,18 @@ void PulseWaveSystem::MergingRiemann(Array<OneD, NekDouble> &Au,
/*
* We solve the six constraint equations via a multivariate Newton
* iteration. Equations are:
- * 1. Backward characteristic: W2(A_R, U_R) = W1(Au_R, Uu_R)
+ * 1. Backward characteristic: W2(A_R, U_R) = W1(Au_R,
+ * Uu_R)
* 2. Forward characteristic 1: W1(A_L1, U_L1) = W1(Au_L1,
* Uu_R1)
* 3. Forward characteristic 2: W1(A_L2, U_L2) = W1(Au_L2,
* Uu_L2)
- * 4. Conservation of mass: Au_R * Uu_R = Au_L1 * Uu_L1 +
- * Au_L2 * Uu_L2
- * 5. Continuity of total pressure 1: rho * Uu_R * Uu_R / 2 + p(Au_R)
- * = rho * Uu_L1 * Uu_L1 / 2 + p(Au_L1)
- * 6. Continuity of total pressure 2: rho * Uu_R * Uu_R / 2 + p(Au_R)
- * = rho * Uu_L2 * Uu_L2 / 2 + p(Au_L2)
+ * 4. Conservation of mass: Au_R * Uu_R = Au_L1 *
+ * Uu_L1 + Au_L2 * Uu_L2
+ * 5. Continuity of total pressure 1: rho * Uu_R * Uu_R / 2 +
+ * p(Au_R) = rho * Uu_L1 * Uu_L1 / 2 + p(Au_L1)
+ * 6. Continuity of total pressure 2: rho * Uu_R * Uu_R / 2 +
+ * p(Au_R) = rho * Uu_L2 * Uu_L2 / 2 + p(Au_L2)
*/
for (size_t i = 0; i < 3; ++i)
{
@@ -1770,8 +1855,8 @@ void PulseWaveSystem::InterfaceRiemann(Array<OneD, NekDouble> &Au,
* 1. Forward characteristic: W1(A_L, U_L) = W1(Au_L, Uu_L)
* 2. Backward characteristic: W2(A_R, U_R) = W2(Au_R, Uu_R)
* 3. Conservation of mass: Au_L * Uu_L = Au_R * Uu_R
- * 4. Continuity of total pressure: rho * Uu_L * Uu_L / 2 + p(Au_L) =
- * rho * Uu_R * Uu_R / 2 + p(Au_R)
+ * 4. Continuity of total pressure: rho * Uu_L * Uu_L / 2 + p(Au_L)
+ * = rho * Uu_R * Uu_R / 2 + p(Au_R)
*/
for (size_t i = 0; i < 2; ++i)
{
@@ -1825,8 +1910,8 @@ void PulseWaveSystem::InterfaceRiemann(Array<OneD, NekDouble> &Au,
void PulseWaveSystem::v_Output(void)
{
/**
- * Write the field data to file. The file is named according to the session
- * name with the extension .fld appended.
+ * Write the field data to file. The file is named according to the
+ * session name with the extension .fld appended.
*/
std::string outname = m_sessionName + ".fld";
@@ -1957,8 +2042,8 @@ NekDouble PulseWaveSystem::v_L2Error(
{
size_t vesselid = field + omega * m_nVariables;
- // since exactsoln is passed for just the first field size we need
- // to reset it each domain
+ // since exactsoln is passed for just the first field size we
+ // need to reset it each domain
Array<OneD, NekDouble> exactsoln(m_vessels[vesselid]->GetNpoints());
m_fields[field] = m_vessels[omega * m_nVariables];
EvaluateExactSolution(field, exactsoln, m_time);
diff --git a/solvers/PulseWaveSolver/EquationSystems/PulseWaveSystem.h b/solvers/PulseWaveSolver/EquationSystems/PulseWaveSystem.h
index 0d3774ec670a2ae5a134cba1ed5791bcdc9ad747..0a2c59ca710ad356396ddbaef47cba475ac04067 100644
--- a/solvers/PulseWaveSolver/EquationSystems/PulseWaveSystem.h
+++ b/solvers/PulseWaveSolver/EquationSystems/PulseWaveSystem.h
@@ -117,6 +117,8 @@ protected:
// keep local copy so can be reordered in parallle
std::map<int, SpatialDomains::CompositeMap> m_domain;
std::vector<int> m_domOrder;
+ std::map<int, std::vector<int>> m_domainToFilterIDs;
+ std::map<int, int> m_filterToVesselID;
Array<OneD, Array<OneD, NekDouble>> m_pressure;
PulseWavePressureAreaSharedPtr m_pressureArea;
diff --git a/solvers/PulseWaveSolver/PulseWaveSolver.cpp b/solvers/PulseWaveSolver/PulseWaveSolver.cpp
index d0f0b2e17380bbd4ec5db7c4b225428596fd7435..79d5370dfc31f1a513c0ac0be771448b4d1eeb5c 100644
--- a/solvers/PulseWaveSolver/PulseWaveSolver.cpp
+++ b/solvers/PulseWaveSolver/PulseWaveSolver.cpp
@@ -64,9 +64,14 @@ int main(int argc, char *argv[])
// Execute driver
drv->Execute();
+ // Print out timings
+ int iolevel = 0;
+
+ session->LoadParameter("IO_Timer_Level", iolevel, -1);
+
// Print out timings
LibUtilities::Timer::PrintElapsedRegions(
- session->GetComm()->GetSpaceComm());
+ session->GetComm()->GetSpaceComm(), std::cout, iolevel);
// Finalise session
session->Finalise();
}
diff --git a/solvers/PulseWaveSolver/Tests/TwoBifurcations.xml b/solvers/PulseWaveSolver/Tests/TwoBifurcations.xml
index 43814f6d2daed356ffecc0946c333eecf2b9511a..b566044a4040449e900a4410dd6ddf23f7f46c43 100644
--- a/solvers/PulseWaveSolver/Tests/TwoBifurcations.xml
+++ b/solvers/PulseWaveSolver/Tests/TwoBifurcations.xml
@@ -108,6 +108,22 @@
<E COMPOSITE="C[8]" NUMMODES="7" FIELDS="A,u" TYPE="MODIFIED" />
</EXPANSIONS>
+ <FILTERS>
+ <FILTER TYPE="HistoryPoints" DOMAIN="0">
+ <PARAM NAME="OutputFile">TimeValues_D0</PARAM>
+ <PARAM NAME="Points">
+ -1.000e+02 0.000e+01 0.000e+00
+ </PARAM>
+ </FILTER>
+ <FILTER TYPE="HistoryPoints" DOMAIN="3">
+ <PARAM NAME="OutputFile">TimeValues_D3</PARAM>
+ <PARAM NAME="Points">
+ 13.000e+01 1.5000e+01 0.000e+00
+ 14.000e+01 1.5000e+01 0.000e+00
+ </PARAM>
+ </FILTER>
+ </FILTERS>
+
<CONDITIONS>
<PARAMETERS>
diff --git a/solvers/PulseWaveSolver/Tests/TwoBifurcations_par.tst b/solvers/PulseWaveSolver/Tests/TwoBifurcations_par.tst
index 9d5fd40e218ae9e13aa308db4aa0b5ecff1eb593..cba600c25cb42933e863b0113d7733461585ff70 100644
--- a/solvers/PulseWaveSolver/Tests/TwoBifurcations_par.tst
+++ b/solvers/PulseWaveSolver/Tests/TwoBifurcations_par.tst
@@ -2,7 +2,7 @@
<test>
<description>Double Bifurcation, P=5</description>
<executable>PulseWaveSolver</executable>
- <parameters>TwoBifurcations.xml</parameters>
+ <parameters>-f TwoBifurcations.xml</parameters>
<processes>2</processes>
<files>
<file description="Session File">TwoBifurcations.xml</file>