Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
FilterBenchmark.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File FilterBenchmark.cpp
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Outputs times when solution crosses a threshold value.
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
37 
38 namespace Nektar
39 {
40 std::string FilterBenchmark::className =
42  "Benchmark",
44 
45 /**
46  * @class FilterBenchmark
47  *
48  * This class records the sequence of activation and repolarisation times across
49  * the entire domain into a two-dimensional storage structure. At each
50  * timestep, the voltage at each point in the domain is examined to identify if
51  * it has crossed the threshold value. If so, the time of crossing is recorded.
52  * Auxiliary arrays hold the current index of each point (i.e. the number of
53  * crossings of the threshold) and the type of the last crossing (activation or
54  * repolarisation).
55  */
56 
57 /**
58  * @param pSession Session reader for IO
59  * @param pParams Parameters of filter
60  */
63  const std::map<std::string, std::string> &pParams)
64  : Filter(pSession)
65 {
66  ASSERTL0(pParams.find("ThresholdValue") != pParams.end(),
67  "Missing parameter 'ThresholdValue'.");
68  m_thresholdValue = atof(pParams.find("ThresholdValue")->second.c_str());
69  ASSERTL0(pParams.find("InitialValue") != pParams.end(),
70  "Missing parameter 'InitialValue'.");
71  m_initialValue = atof(pParams.find("InitialValue")->second.c_str());
72  ASSERTL0(!(pParams.find("OutputFile")->second.empty()),
73  "Missing parameter 'OutputFile'.");
74  m_outputFile = pParams.find("OutputFile")->second;
75 
76  m_startTime = 0.0;
77  if (pParams.find("StartTime") != pParams.end())
78  {
79  m_startTime = atof(pParams.find("StartTime")->second.c_str());
80  }
81 
83  pSession->GetComm());
84 
85 }
86 
87 
88 /**
89  *
90  */
92 {
93 
94 }
95 
96 
97 /*
98  * Initialises the storage.
99  * @param pFields Field storage expansion lists
100  * @param time Current time
101  */
103  const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields,
104  const NekDouble &time)
105 {
106  m_threshold.push_back(Array<OneD, NekDouble>(
107  pFields[0]->GetNpoints(), m_initialValue));
108 
109  m_idx = Array<OneD, int> (pFields[0]->GetNpoints(), 0);
110  m_polarity = Array<OneD, int> (pFields[0]->GetNpoints(), -1);
111 }
112 
113 
114 /**
115  * Checks each point in the domain to determine if it has crossed the threshold.
116  * The direction of crossing is determined. Additional storage is allocated if
117  * needed.
118  * @param pFields Field storage expansion lists
119  * @param time Current time
120  */
122  const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields,
123  const NekDouble &time)
124 {
125  // Only proceed if the start time has passed
126  if (time < m_startTime)
127  {
128  return;
129  }
130 
131  // Examine each point in turn
132  for (int i = 0; i < pFields[0]->GetNpoints(); ++i)
133  {
134  if ((m_polarity[i] == -1 &&
135  pFields[0]->GetPhys()[i] > m_thresholdValue) ||
136  (m_polarity[i] == 1 &&
137  pFields[0]->GetPhys()[i] < m_thresholdValue))
138  {
139  // If APD less than 50ms, remove last activation
140  if (m_polarity[i] == 1 &&
141  time - m_threshold[m_idx[i]][i] < 50)
142  {
143  m_idx[i]--;
145  }
146  else
147  {
148  m_threshold[m_idx[i]][i] = time;
149  m_idx[i]++;
150  }
151  // Update polarity of last crossing
152  m_polarity[i] *= -1;
153  }
154  }
155 
156  // Allocate additional storage if any point has as many crossings as
157  // current storage permits.
158  int max_idx = Vmath::Vmax(pFields[0]->GetNpoints(), m_idx, 1);
159  pFields[0]->GetSession()->GetComm()->AllReduce(max_idx,
161  if (m_threshold.size() == max_idx)
162  {
163  m_threshold.push_back(Array<OneD, NekDouble>(
164  pFields[0]->GetNpoints(), m_initialValue));
165  }
166 
167 }
168 
169 
170 /**
171  * Writes out the crossings to file.
172  * @param pFields Field storage expansion list.
173  * @param time Current time.
174  */
176  const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields,
177  const NekDouble &time)
178 {
179  for (int i = 0; i < m_threshold.size() - 1; ++i)
180  {
181  std::stringstream vOutputFilename;
182  vOutputFilename << m_outputFile << "_" << i << ".fld";
183 
184  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef
185  = pFields[0]->GetFieldDefinitions();
186  std::vector<std::vector<NekDouble> > FieldData(FieldDef.size());
187 
188  Array<OneD, NekDouble> vCoeffs(pFields[0]->GetNcoeffs());
189  pFields[0]->FwdTrans_IterPerExp(m_threshold[i], vCoeffs);
190 
191  // copy Data into FieldData and set variable
192  for(int i = 0; i < FieldDef.size(); ++i)
193  {
194  // Could do a search here to find correct variable
195  FieldDef[i]->m_fields.push_back("m");
196  pFields[0]->AppendFieldData(FieldDef[i], FieldData[i], vCoeffs);
197  }
198 
199  m_fld->Write(vOutputFilename.str(),FieldDef,FieldData);
200  }
201 }
202 
203 
204 /**
205  * @return This filter is time dependent.
206  */
208 {
209  return true;
210 }
211 
212 }