Nektar++
FilterMovingBody.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File FilterMovingBody.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: Output moving body forces during time-stepping.
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <LibUtilities/Memory/NekMemoryManager.hpp>
37 #include <iomanip>
38 #include <LocalRegions/Expansion1D.h>
39 #include <LocalRegions/Expansion2D.h>
40 #include <LocalRegions/Expansion3D.h>
41 #include <IncNavierStokesSolver/Filters/FilterMovingBody.h>
42 
43 namespace Nektar
44 {
45 
46 std::string FilterMovingBody::className = SolverUtils::GetFilterFactory().
47  RegisterCreatorFunction("MovingBody",
48  FilterMovingBody::create,
49  "Moving Body Filter");
50 /**
51  *
52  */
53 FilterMovingBody::FilterMovingBody(
54  const LibUtilities::SessionReaderSharedPtr &pSession,
55  const std::map<std::string, std::string> &pParams)
56  : Filter(pSession),
57  m_session(pSession)
58 {
59  if (pParams.find("OutputFile") == pParams.end())
60  {
61  m_outputFile_fce = pSession->GetSessionName();
62  m_outputFile_mot = pSession->GetSessionName();
63  }
64  else
65  {
66  ASSERTL0(!(pParams.find("OutputFile")->second.empty()),
67  "Missing parameter 'OutputFile'.");
68 
69  m_outputFile_fce = pParams.find("OutputFile")->second;
70  m_outputFile_mot = pParams.find("OutputFile")->second;
71  }
72  if (!(m_outputFile_fce.length() >= 4 &&
73  m_outputFile_fce.substr(m_outputFile_fce.length() - 4) == ".fce"))
74  {
75  m_outputFile_fce += ".fce";
76  }
77 
78  if (!(m_outputFile_mot.length() >= 4 &&
79  m_outputFile_mot.substr(m_outputFile_mot.length() - 4) == ".mot"))
80  {
81  m_outputFile_mot += ".mot";
82  }
83  if (pParams.find("OutputFrequency") == pParams.end())
84  {
85  m_outputFrequency = 1;
86  }
87  else
88  {
89  m_outputFrequency = atoi(
90  pParams.find("OutputFrequency")->second.c_str());
91  }
92 
93  pSession->MatchSolverInfo("Homogeneous", "1D", m_isHomogeneous1D, false);
94 
95  ASSERTL0(m_isHomogeneous1D, "Moving Body implemented just for 3D "
96  "Homogeneous 1D discetisations.");
97 
98  //specify the boundary to calculate the forces
99  if (pParams.find("Boundary") == pParams.end())
100  {
101  ASSERTL0(false, "Missing parameter 'Boundary'.");
102  }
103  else
104  {
105  ASSERTL0(!(pParams.find("Boundary")->second.empty()),
106  "Missing parameter 'Boundary'.");
107  m_BoundaryString = pParams.find("Boundary")->second;
108  }
109 }
110 
111 
112 /**
113  *
114  */
115 FilterMovingBody::~FilterMovingBody()
116 {
117 
118 }
119 
120 
121 /**
122  *
123  */
124 void FilterMovingBody::v_Initialise(
125  const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields,
126  const NekDouble &time)
127 {
128  m_index = 0;
129  m_outputStream = Array<OneD, std::ofstream>(2);
130  // Parse the boundary regions into a list.
131  std::string::size_type FirstInd = m_BoundaryString.find_first_of('[') + 1;
132  std::string::size_type LastInd = m_BoundaryString.find_last_of(']') - 1;
133 
134  ASSERTL0(FirstInd <= LastInd,
135  (std::string("Error reading boundary region definition:") +
136  m_BoundaryString).c_str());
137 
138  std::string IndString = m_BoundaryString.substr(FirstInd,
139  LastInd - FirstInd + 1);
140 
141  bool parseGood = ParseUtils::GenerateSeqVector(IndString.c_str(),
142  m_boundaryRegionsIdList);
143 
144  ASSERTL0(parseGood && !m_boundaryRegionsIdList.empty(),
145  (std::string("Unable to read boundary regions index "
146  "range for FilterAeroForces: ") + IndString).c_str());
147 
148  // determine what boundary regions need to be considered
149  int cnt;
150 
151  unsigned int numBoundaryRegions
152  = pFields[0]->GetBndConditions().num_elements();
153 
154  m_boundaryRegionIsInList.insert(m_boundaryRegionIsInList.end(),
155  numBoundaryRegions, 0);
156 
157  SpatialDomains::BoundaryConditions bcs(m_session,pFields[0]->GetGraph());
158 
159  const SpatialDomains::BoundaryRegionCollection &bregions
160  = bcs.GetBoundaryRegions();
161 
162  SpatialDomains::BoundaryRegionCollection::const_iterator it;
163 
164  for (cnt = 0, it = bregions.begin(); it != bregions.end();
165  ++it, cnt++)
166  {
167  if ( std::find(m_boundaryRegionsIdList.begin(),
168  m_boundaryRegionsIdList.end(), it->first) !=
169  m_boundaryRegionsIdList.end() )
170  {
171  m_boundaryRegionIsInList[cnt] = 1;
172  }
173  }
174 
175  LibUtilities::CommSharedPtr vComm = pFields[0]->GetComm();
176 
177  if (vComm->GetRank() == 0)
178  {
179  // Open output stream for cable forces
180  m_outputStream[0].open(m_outputFile_fce.c_str());
181  m_outputStream[0] << "#";
182  m_outputStream[0].width(7);
183  m_outputStream[0] << "Time";
184  m_outputStream[0].width(25);
185  m_outputStream[0] << "z";
186  m_outputStream[0].width(25);
187  m_outputStream[0] << "Fx (press)";
188  m_outputStream[0].width(25);
189  m_outputStream[0] << "Fx (visc)";
190  m_outputStream[0].width(25);
191  m_outputStream[0] << "Fx (tot)";
192  m_outputStream[0].width(25);
193  m_outputStream[0] << "Fy (press)";
194  m_outputStream[0].width(25);
195  m_outputStream[0] << "Fy (visc)";
196  m_outputStream[0].width(25);
197  m_outputStream[0] << "Fy (tot)";
198  m_outputStream[0] << endl;
199 
200  // Open output stream for cable motions
201  m_outputStream[1].open(m_outputFile_mot.c_str());
202  m_outputStream[1] << "#";
203  m_outputStream[1].width(7);
204  m_outputStream[1] << "Time";
205  m_outputStream[1].width(25);
206  m_outputStream[1] << "z";
207  m_outputStream[1].width(25);
208  m_outputStream[1] << "Disp_x";
209  m_outputStream[1].width(25);
210  m_outputStream[1] << "Vel_x";
211  m_outputStream[1].width(25);
212  m_outputStream[1] << "Acel_x";
213  m_outputStream[1].width(25);
214  m_outputStream[1] << "Disp_y";
215  m_outputStream[1].width(25);
216  m_outputStream[1] << "Vel_y";
217  m_outputStream[1].width(25);
218  m_outputStream[1] << "Acel_y";
219  m_outputStream[1] << endl;
220  }
221 }
222 
223 
224 /**
225  *
226  */
227 void FilterMovingBody::UpdateForce(
228  const LibUtilities::SessionReaderSharedPtr &pSession,
229  const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields,
230  Array<OneD, NekDouble> &Aeroforces,
231  const NekDouble &time)
232 {
233  int n, cnt, elmtid, nq, offset, boundary;
234  int nt = pFields[0]->GetNpoints();
235  int dim = pFields.num_elements()-1;
236 
237  StdRegions::StdExpansionSharedPtr elmt;
238  Array<OneD, int> BoundarytoElmtID;
239  Array<OneD, int> BoundarytoTraceID;
240  Array<OneD, MultiRegions::ExpListSharedPtr> BndExp;
241 
242  Array<OneD, const NekDouble> P(nt);
243  Array<OneD, const NekDouble> U(nt);
244  Array<OneD, const NekDouble> V(nt);
245  Array<OneD, const NekDouble> W(nt);
246 
247  Array<OneD, Array<OneD, NekDouble> > gradU(dim);
248  Array<OneD, Array<OneD, NekDouble> > gradV(dim);
249  Array<OneD, Array<OneD, NekDouble> > gradW(dim);
250 
251  Array<OneD, Array<OneD, NekDouble> > fgradU(dim);
252  Array<OneD, Array<OneD, NekDouble> > fgradV(dim);
253  Array<OneD, Array<OneD, NekDouble> > fgradW(dim);
254 
255  LibUtilities::CommSharedPtr vComm = pFields[0]->GetComm();
256  LibUtilities::CommSharedPtr vRowComm = vComm->GetRowComm();
257  LibUtilities::CommSharedPtr vColComm = vComm->GetColumnComm();
258  LibUtilities::CommSharedPtr vCColComm = vColComm->GetColumnComm();
259 
260  // set up storage space for forces on all the planes (z-positions)
261  // on each processors some of them will remain empty as we may
262  // have just few planes per processor
263  int Num_z_pos = pFields[0]->GetHomogeneousBasis()->GetNumModes();
264  Array<OneD, NekDouble> Fx(Num_z_pos,0.0);
265  Array<OneD, NekDouble> Fxp(Num_z_pos,0.0);
266  Array<OneD, NekDouble> Fxv(Num_z_pos,0.0);
267  Array<OneD, NekDouble> Fy(Num_z_pos,0.0);
268  Array<OneD, NekDouble> Fyp(Num_z_pos,0.0);
269  Array<OneD, NekDouble> Fyv(Num_z_pos,0.0);
270 
271 
272  NekDouble rho = (pSession->DefinesParameter("rho"))
273  ? (pSession->GetParameter("rho"))
274  : 1;
275  NekDouble mu = rho*pSession->GetParameter("Kinvis");
276 
277  for(int i = 0; i < pFields.num_elements(); ++i)
278  {
279  pFields[i]->SetWaveSpace(false);
280  pFields[i]->BwdTrans(pFields[i]->GetCoeffs(),
281  pFields[i]->UpdatePhys());
282  pFields[i]->SetPhysState(true);
283  }
284 
285  // Get the number of local planes on the process and their IDs
286  // to properly locate the forces in the Fx, Fy etc. vectors.
287  Array<OneD, unsigned int> ZIDs;
288  ZIDs = pFields[0]->GetZIDs();
289  int local_planes = ZIDs.num_elements();
290 
291  // Homogeneous 1D case Compute forces on all WALL boundaries
292  // This only has to be done on the zero (mean) Fourier mode.
293  for(int plane = 0 ; plane < local_planes; plane++)
294  {
295  pFields[0]->GetPlane(plane)->GetBoundaryToElmtMap(BoundarytoElmtID,
296  BoundarytoTraceID);
297  BndExp = pFields[0]->GetPlane(plane)->GetBndCondExpansions();
298  StdRegions::StdExpansionSharedPtr bc;
299 
300  // loop over the types of boundary conditions
301  for(cnt = n = 0; n < BndExp.num_elements(); ++n)
302  {
303  if(m_boundaryRegionIsInList[n] == 1)
304  {
305  for(int i = 0; i < BndExp[n]->GetExpSize(); ++i, cnt++)
306  {
307  // find element of this expansion.
308  elmtid = BoundarytoElmtID[cnt];
309  elmt = pFields[0]->GetPlane(plane)->GetExp(elmtid);
310  nq = elmt->GetTotPoints();
311  offset = pFields[0]->GetPlane(plane)
312  ->GetPhys_Offset(elmtid);
313 
314  // Initialise local arrays for the velocity
315  // gradients size of total number of quadrature
316  // points for each element (hence local).
317  for(int j = 0; j < dim; ++j)
318  {
319  gradU[j] = Array<OneD, NekDouble>(nq,0.0);
320  gradV[j] = Array<OneD, NekDouble>(nq,0.0);
321  gradW[j] = Array<OneD, NekDouble>(nq,0.0);
322  }
323 
324  // identify boundary of element
325  boundary = BoundarytoTraceID[cnt];
326 
327  // Extract fields
328  U = pFields[0]->GetPlane(plane)->GetPhys() + offset;
329  V = pFields[1]->GetPlane(plane)->GetPhys() + offset;
330  P = pFields[3]->GetPlane(plane)->GetPhys() + offset;
331 
332  // compute the gradients
333  elmt->PhysDeriv(U,gradU[0],gradU[1]);
334  elmt->PhysDeriv(V,gradV[0],gradV[1]);
335 
336  // Get face 1D expansion from element expansion
337  bc = BndExp[n]->GetExp(i)->as<LocalRegions::Expansion1D>();
338 
339  // number of points on the boundary
340  int nbc = bc->GetTotPoints();
341 
342  // several vectors for computing the forces
343  Array<OneD, NekDouble> Pb(nbc,0.0);
344 
345  for(int j = 0; j < dim; ++j)
346  {
347  fgradU[j] = Array<OneD, NekDouble>(nbc,0.0);
348  fgradV[j] = Array<OneD, NekDouble>(nbc,0.0);
349  }
350 
351  Array<OneD, NekDouble> drag_t(nbc,0.0);
352  Array<OneD, NekDouble> lift_t(nbc,0.0);
353  Array<OneD, NekDouble> drag_p(nbc,0.0);
354  Array<OneD, NekDouble> lift_p(nbc,0.0);
355  Array<OneD, NekDouble> temp(nbc,0.0);
356  Array<OneD, NekDouble> temp2(nbc,0.0);
357 
358  // identify boundary of element .
359  boundary = BoundarytoTraceID[cnt];
360 
361  // extraction of the pressure and wss on the
362  // boundary of the element
363  elmt->GetEdgePhysVals(boundary,bc,P,Pb);
364 
365  for(int j = 0; j < dim; ++j)
366  {
367  elmt->GetEdgePhysVals(boundary,bc,gradU[j],fgradU[j]);
368  elmt->GetEdgePhysVals(boundary,bc,gradV[j],fgradV[j]);
369  }
370 
371  //normals of the element
372  const Array<OneD, Array<OneD, NekDouble> > &normals
373  = elmt->GetEdgeNormal(boundary);
374 
375  //
376  // Compute viscous tractive forces on wall from
377  //
378  // t_i = - T_ij * n_j (minus sign for force
379  // exerted BY fluid ON wall),
380  //
381  // where
382  //
383  // T_ij = viscous stress tensor (here in Cartesian
384  // coords)
385  // dU_i dU_j
386  // = RHO * KINVIS * ( ---- + ---- ) .
387  // dx_j dx_i
388 
389  //a) DRAG TERMS
390  //-rho*kinvis*(2*du/dx*nx+(du/dy+dv/dx)*ny)
391 
392  Vmath::Vadd(nbc, fgradU[1], 1, fgradV[0], 1, drag_t, 1);
393  Vmath::Vmul(nbc, drag_t, 1, normals[1], 1, drag_t, 1);
394 
395  Vmath::Smul(nbc, 2.0, fgradU[0], 1, fgradU[0], 1);
396  Vmath::Vmul(nbc, fgradU[0], 1, normals[0], 1, temp2, 1);
397  Vmath::Smul(nbc, 0.5, fgradU[0], 1, fgradU[0], 1);
398 
399  Vmath::Vadd(nbc, temp2, 1, drag_t, 1, drag_t, 1);
400  Vmath::Smul(nbc, -mu, drag_t, 1, drag_t, 1);
401 
402  //zero temporary storage vector
403  Vmath::Zero(nbc, temp, 0.0);
404  Vmath::Zero(nbc, temp2, 0.0);
405 
406 
407  //b) LIFT TERMS
408  //-rho*kinvis*(2*dv/dy*ny+(du/dy+dv/dx)*nx)
409 
410  Vmath::Vadd(nbc, fgradU[1], 1, fgradV[0], 1, lift_t, 1);
411  Vmath::Vmul(nbc, lift_t, 1, normals[0], 1, lift_t, 1);
412 
413  Vmath::Smul(nbc, 2.0, fgradV[1], 1, fgradV[1], 1);
414  Vmath::Vmul(nbc, fgradV[1], 1, normals[1], 1, temp2, 1);
415  Vmath::Smul(nbc, -0.5, fgradV[1], 1, fgradV[1], 1);
416 
417 
418  Vmath::Vadd(nbc, temp2, 1, lift_t, 1, lift_t, 1);
419  Vmath::Smul(nbc, -mu, lift_t, 1, lift_t, 1);
420 
421  // Compute normal tractive forces on all WALL
422  // boundaries
423 
424  Vmath::Vvtvp(nbc, Pb, 1, normals[0], 1, drag_p, 1,
425  drag_p, 1);
426  Vmath::Vvtvp(nbc, Pb, 1, normals[1], 1, lift_p, 1,
427  lift_p, 1);
428 
429  //integration over the boundary
430  Fxv[ZIDs[plane]] += bc->Integral(drag_t);
431  Fyv[ZIDs[plane]] += bc->Integral(lift_t);
432 
433  Fxp[ZIDs[plane]] += bc->Integral(drag_p);
434  Fyp[ZIDs[plane]] += bc->Integral(lift_p);
435  }
436  }
437  else
438  {
439  cnt += BndExp[n]->GetExpSize();
440  }
441  }
442  }
443 
444  for(int i = 0; i < pFields.num_elements(); ++i)
445  {
446  pFields[i]->SetWaveSpace(true);
447  pFields[i]->BwdTrans(pFields[i]->GetCoeffs(),
448  pFields[i]->UpdatePhys());
449  pFields[i]->SetPhysState(false);
450  }
451 
452  // In case we are using an hybrid parallelisation we need
453  // to reduce the forces on the same plane coming from
454  // different mesh partitions.
455  // It is quite an expensive communication, therefore
456  // we check to make sure it is actually required.
457 
458  if(vComm->GetRowComm()->GetSize() > 0)
459  {
460  // NOTE 1: We can eventually sum the viscous and pressure
461  // component before doing the communication, thus
462  // reducing by a factor 2 the communication.
463  // NOTE 2: We may want to set up in the Comm class an AllReduce
464  // routine wich can handle arrays more efficiently
465  for(int plane = 0 ; plane < local_planes; plane++)
466  {
467  vRowComm->AllReduce(Fxp[ZIDs[plane]], LibUtilities::ReduceSum);
468  vRowComm->AllReduce(Fxv[ZIDs[plane]], LibUtilities::ReduceSum);
469  vRowComm->AllReduce(Fyp[ZIDs[plane]], LibUtilities::ReduceSum);
470  vRowComm->AllReduce(Fyv[ZIDs[plane]], LibUtilities::ReduceSum);
471  }
472  }
473 
474  // At this point on rank (0,n) of the Mesh partion communicator we have
475  // the total areo forces on the planes which are on the same column
476  // communicator. Since the planes are scattered on different processors
477  // some of the entries of the vector Fxp, Fxp etc. are still zero.
478  // Now we need to reduce the values on a single vector on rank (0,0) of the
479  // global communicator.
480  if(!pSession->DefinesSolverInfo("HomoStrip"))
481  {
482  if(vComm->GetRowComm()->GetRank() == 0)
483  {
484  for(int z = 0 ; z < Num_z_pos; z++)
485  {
486  vColComm->AllReduce(Fxp[z], LibUtilities::ReduceSum);
487  vColComm->AllReduce(Fxv[z], LibUtilities::ReduceSum);
488  vColComm->AllReduce(Fyp[z], LibUtilities::ReduceSum);
489  vColComm->AllReduce(Fyv[z], LibUtilities::ReduceSum);
490  }
491  }
492  }
493  else
494  {
495  if(vComm->GetRowComm()->GetRank() == 0)
496  {
497  for(int z = 0 ; z < Num_z_pos; z++)
498  {
499  vCColComm->AllReduce(Fxp[z], LibUtilities::ReduceSum);
500  vCColComm->AllReduce(Fxv[z], LibUtilities::ReduceSum);
501  vCColComm->AllReduce(Fyp[z], LibUtilities::ReduceSum);
502  vCColComm->AllReduce(Fyv[z], LibUtilities::ReduceSum);
503  }
504  }
505  }
506 
507  if(!pSession->DefinesSolverInfo("HomoStrip"))
508  {
509  //set the forces imparted on the cable's wall
510  for(int plane = 0 ; plane < local_planes; plane++)
511  {
512  Aeroforces[plane] = Fxp[ZIDs[plane]]
513  + Fxv[ZIDs[plane]];
514  Aeroforces[plane + local_planes] = Fyp[ZIDs[plane]]
515  + Fyv[ZIDs[plane]];
516  }
517 
518  // Only output every m_outputFrequency.
519  if ((m_index++) % m_outputFrequency)
520  {
521  return;
522  }
523 
524  // At thi point in rank (0,0) we have the full vectors
525  // containing Fxp,Fxv,Fyp and Fyv where different positions
526  // in the vectors correspond to different planes.
527  // Here we write it to file. We do it just on one porcess
528 
529  Array<OneD, NekDouble> z_coords(Num_z_pos,0.0);
530  Array<OneD, const NekDouble> pts
531  = pFields[0]->GetHomogeneousBasis()->GetZ();
532 
533  NekDouble LZ;
534  pSession->LoadParameter("LZ", LZ);
535  Vmath::Smul(Num_z_pos,LZ/2.0,pts,1,z_coords,1);
536  Vmath::Sadd(Num_z_pos,LZ/2.0,z_coords,1,z_coords,1);
537  if (vComm->GetRank() == 0)
538  {
539 
540  Vmath::Vadd(Num_z_pos,Fxp,1,Fxv,1,Fx,1);
541  Vmath::Vadd(Num_z_pos,Fyp,1,Fyv,1,Fy,1);
542 
543  for(int i = 0 ; i < Num_z_pos; i++)
544  {
545  m_outputStream[0].width(8);
546  m_outputStream[0] << setprecision(6) << time;
547 
548  m_outputStream[0].width(25);
549  m_outputStream[0] << setprecision(6) << z_coords[i];
550 
551  m_outputStream[0].width(25);
552  m_outputStream[0] << setprecision(8) << Fxp[i];
553  m_outputStream[0].width(25);
554  m_outputStream[0] << setprecision(8) << Fxv[i];
555  m_outputStream[0].width(25);
556  m_outputStream[0] << setprecision(8) << Fx[i];
557 
558  m_outputStream[0].width(25);
559  m_outputStream[0] << setprecision(8) << Fyp[i];
560  m_outputStream[0].width(25);
561  m_outputStream[0] << setprecision(8) << Fyv[i];
562  m_outputStream[0].width(25);
563  m_outputStream[0] << setprecision(8) << Fy[i];
564  m_outputStream[0] << endl;
565  }
566  }
567  }
568  else
569  {
570  // average the forces over local planes for each processor
571  Array<OneD, NekDouble> fces(6,0.0);
572  for(int plane = 0 ; plane < local_planes; plane++)
573  {
574  fces[0] += Fxp[ZIDs[plane]] + Fxv[ZIDs[plane]];
575  fces[1] += Fyp[ZIDs[plane]] + Fyv[ZIDs[plane]];
576  fces[2] += Fxp[ZIDs[plane]] ;
577  fces[3] += Fyp[ZIDs[plane]] ;
578  fces[4] += Fxv[ZIDs[plane]] ;
579  fces[5] += Fyv[ZIDs[plane]] ;
580  }
581 
582  fces[0] = fces[0]/local_planes;
583  fces[1] = fces[1]/local_planes;
584  fces[2] = fces[2]/local_planes;
585  fces[3] = fces[3]/local_planes;
586  fces[4] = fces[4]/local_planes;
587  fces[5] = fces[5]/local_planes;
588 
589  // average the forces over communicators within each strip
590  vCColComm->AllReduce(fces[0], LibUtilities::ReduceSum);
591  vCColComm->AllReduce(fces[1], LibUtilities::ReduceSum);
592  vCColComm->AllReduce(fces[2], LibUtilities::ReduceSum);
593  vCColComm->AllReduce(fces[3], LibUtilities::ReduceSum);
594  vCColComm->AllReduce(fces[4], LibUtilities::ReduceSum);
595  vCColComm->AllReduce(fces[5], LibUtilities::ReduceSum);
596 
597  int npts = vComm->GetColumnComm()->GetColumnComm()->GetSize();
598 
599  fces[0] = fces[0]/npts;
600  fces[1] = fces[1]/npts;
601  fces[2] = fces[2]/npts;
602  fces[3] = fces[3]/npts;
603  fces[4] = fces[4]/npts;
604  fces[5] = fces[5]/npts;
605 
606  for(int plane = 0 ; plane < local_planes; plane++)
607  {
608  Aeroforces[plane] = fces[0];
609  Aeroforces[plane+local_planes] = fces[1];
610  }
611 
612  // Only output every m_outputFrequency.
613  if ((m_index) % m_outputFrequency)
614  {
615  return;
616  }
617 
618  int colrank = vColComm->GetRank();
619  int nstrips;
620 
621  NekDouble DistStrip;
622 
623  pSession->LoadParameter("Strip_Z", nstrips);
624  pSession->LoadParameter("DistStrip", DistStrip);
625 
626  Array<OneD, NekDouble> z_coords(nstrips);
627  for(int i = 0; i < nstrips; i++)
628  {
629  z_coords[i] = i * DistStrip;
630  }
631 
632  if(colrank == 0)
633  {
634  m_outputStream[0].width(8);
635  m_outputStream[0] << setprecision(6) << time;
636 
637  m_outputStream[0].width(25);
638  m_outputStream[0] << setprecision(6) << z_coords[0];
639 
640  m_outputStream[0].width(25);
641  m_outputStream[0] << setprecision(8) << fces[2];
642  m_outputStream[0].width(25);
643  m_outputStream[0] << setprecision(8) << fces[4];
644  m_outputStream[0].width(25);
645  m_outputStream[0] << setprecision(8) << fces[0];
646 
647  m_outputStream[0].width(25);
648  m_outputStream[0] << setprecision(8) << fces[3];
649  m_outputStream[0].width(25);
650  m_outputStream[0] << setprecision(8) << fces[5];
651  m_outputStream[0].width(25);
652  m_outputStream[0] << setprecision(8) << fces[1];
653  m_outputStream[0] << endl;
654 
655  for(int i = 1; i < nstrips; i++)
656  {
657  vColComm->Recv(i, fces);
658 
659  m_outputStream[0].width(8);
660  m_outputStream[0] << setprecision(6) << time;
661 
662  m_outputStream[0].width(25);
663  m_outputStream[0] << setprecision(6) << z_coords[i];
664 
665  m_outputStream[0].width(25);
666  m_outputStream[0] << setprecision(8) << fces[2];
667  m_outputStream[0].width(25);
668  m_outputStream[0] << setprecision(8) << fces[4];
669  m_outputStream[0].width(25);
670  m_outputStream[0] << setprecision(8) << fces[0];
671 
672  m_outputStream[0].width(25);
673  m_outputStream[0] << setprecision(8) << fces[3];
674  m_outputStream[0].width(25);
675  m_outputStream[0] << setprecision(8) << fces[5];
676  m_outputStream[0].width(25);
677  m_outputStream[0] << setprecision(8) << fces[1];
678  m_outputStream[0] << endl;
679  }
680  }
681  else
682  {
683  for(int i = 1; i < nstrips; i++)
684  {
685  if(colrank == i)
686  {
687  vColComm->Send(0, fces);
688  }
689  }
690  }
691  }
692 }
693 
694 
695 /**
696  *
697  */
698 void FilterMovingBody::UpdateMotion(
699  const LibUtilities::SessionReaderSharedPtr &pSession,
700  const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields,
701  Array<OneD, NekDouble> &MotionVars,
702  const NekDouble &time)
703 {
704  // Only output every m_outputFrequency.
705  if ((m_index++) % m_outputFrequency)
706  {
707  return;
708  }
709 
710  // Get the number of local planes on the process and their IDs
711  // to properly locate the forces in the Fx, Fy etc. vectors.
712  Array<OneD, unsigned int> ZIDs;
713  ZIDs = pFields[0]->GetZIDs();
714  int local_planes = ZIDs.num_elements();
715 
716  LibUtilities::CommSharedPtr vColComm
717  = pFields[0]->GetComm()->GetColumnComm();
718 
719  //
720  if(!pSession->DefinesSolverInfo("HomoStrip"))
721  {
722  int Num_z_pos = pFields[0]->GetHomogeneousBasis()->GetNumModes();
723  Array<OneD, NekDouble> z_coords(Num_z_pos,0.0);
724  Array<OneD, const NekDouble> pts
725  = pFields[0]->GetHomogeneousBasis()->GetZ();
726 
727  NekDouble LZ;
728  pSession->LoadParameter("LZ", LZ);
729  Vmath::Smul(Num_z_pos,LZ/2.0,pts,1,z_coords,1);
730  Vmath::Sadd(Num_z_pos,LZ/2.0,z_coords,1,z_coords,1);
731 
732  //get and output moving body variables
733  int nStrVars = 3;
734  Array<OneD, Array<OneD, NekDouble> > Motion_x(nStrVars);
735  Array<OneD, Array<OneD, NekDouble> > Motion_y(nStrVars);
736 
737  for(int i = 0; i < nStrVars; i++)
738  {
739  Motion_x[i] = Array<OneD, NekDouble>(local_planes,0.0);
740  Motion_y[i] = Array<OneD, NekDouble>(local_planes,0.0);
741  }
742 
743  for(int plane = 0; plane < local_planes; plane++)
744  {
745  for (int var = 0; var < nStrVars; var++)
746  {
747  int xoffset = var*local_planes+plane;
748  int yoffset = nStrVars*local_planes+xoffset;
749  Motion_x[var][plane] = MotionVars[xoffset];
750  Motion_y[var][plane] = MotionVars[yoffset];
751  }
752  }
753 
754  Array <OneD, NekDouble> CableAccelX;
755  Array <OneD, NekDouble> CableVelocX;
756  Array <OneD, NekDouble> CableDisplX;
757  Array <OneD, NekDouble> CableAccelY;
758  Array <OneD, NekDouble> CableVelocY;
759  Array <OneD, NekDouble> CableDisplY;
760 
761  int npoints = Motion_x[0].num_elements();
762  CableAccelX = Array <OneD, NekDouble>(npoints);
763  CableVelocX = Array <OneD, NekDouble>(npoints);
764  CableDisplX = Array <OneD, NekDouble>(npoints);
765  CableAccelY = Array <OneD, NekDouble>(npoints);
766  CableVelocY = Array <OneD, NekDouble>(npoints);
767  CableDisplY = Array <OneD, NekDouble>(npoints);
768 
769  Vmath::Vcopy(npoints, Motion_x[0], 1, CableDisplX, 1);
770  Vmath::Vcopy(npoints, Motion_x[1], 1, CableVelocX, 1);
771  Vmath::Vcopy(npoints, Motion_x[2], 1, CableAccelX, 1);
772  Vmath::Vcopy(npoints, Motion_y[0], 1, CableDisplY, 1);
773  Vmath::Vcopy(npoints, Motion_y[1], 1, CableVelocY, 1);
774  Vmath::Vcopy(npoints, Motion_y[2], 1, CableAccelY, 1);
775 
776  int colrank = vColComm->GetRank();
777  int nproc = vColComm->GetSize();
778  // Send to root process.
779  if (colrank == 0)
780  {
781  for (int j = 0; j <Motion_x[0].num_elements(); j++)
782  {
783  m_outputStream[1].width(8);
784  m_outputStream[1] << setprecision(6) << time;
785  m_outputStream[1].width(25);
786  m_outputStream[1] << setprecision(6) << z_coords[j];
787  m_outputStream[1].width(25);
788  m_outputStream[1] << setprecision(8) << CableDisplX[j];
789  m_outputStream[1].width(25);
790  m_outputStream[1] << setprecision(8) << CableVelocX[j];
791  m_outputStream[1].width(25);
792  m_outputStream[1] << setprecision(8) << CableAccelX[j];
793  m_outputStream[1].width(25);
794  m_outputStream[1] << setprecision(8) << CableDisplY[j];
795  m_outputStream[1].width(25);
796  m_outputStream[1] << setprecision(8) << CableVelocY[j];
797  m_outputStream[1].width(25);
798  m_outputStream[1] << setprecision(8) << CableAccelY[j];
799  m_outputStream[1] << endl;
800  }
801 
802  for (int i = 1; i < nproc; ++i)
803  {
804  vColComm->Recv(i, CableAccelX);
805  vColComm->Recv(i, CableVelocX);
806  vColComm->Recv(i, CableDisplX);
807  vColComm->Recv(i, CableAccelY);
808  vColComm->Recv(i, CableVelocY);
809  vColComm->Recv(i, CableDisplY);
810 
811 
812  for (int j = 0; j < Motion_x[0].num_elements(); ++j)
813  {
814  int n = Num_z_pos/nproc * i + j;
815  m_outputStream[1].width(8);
816  m_outputStream[1] << setprecision(6) << time;
817  m_outputStream[1].width(25);
818  m_outputStream[1] << setprecision(6) << z_coords[n];
819  m_outputStream[1].width(25);
820  m_outputStream[1] << setprecision(8) << CableDisplX[j];
821  m_outputStream[1].width(25);
822  m_outputStream[1] << setprecision(8) << CableVelocX[j];
823  m_outputStream[1].width(25);
824  m_outputStream[1] << setprecision(8) << CableAccelX[j];
825  m_outputStream[1].width(25);
826  m_outputStream[1] << setprecision(8) << CableDisplY[j];
827  m_outputStream[1].width(25);
828  m_outputStream[1] << setprecision(8) << CableVelocY[j];
829  m_outputStream[1].width(25);
830  m_outputStream[1] << setprecision(8) << CableAccelY[j];
831  m_outputStream[1] << endl;
832  }
833  }
834  }
835  else
836  {
837  vColComm->Send(0, CableAccelX);
838  vColComm->Send(0, CableVelocX);
839  vColComm->Send(0, CableDisplX);
840  vColComm->Send(0, CableAccelY);
841  vColComm->Send(0, CableVelocY);
842  vColComm->Send(0, CableDisplY);
843  }
844  }
845  else
846  {
847  int colrank = vColComm->GetRank();
848  int nstrips;
849 
850  NekDouble DistStrip;
851 
852  pSession->LoadParameter("Strip_Z", nstrips);
853  pSession->LoadParameter("DistStrip", DistStrip);
854 
855  Array<OneD, NekDouble> z_coords(nstrips);
856  for(int i = 0; i < nstrips; i++)
857  {
858  z_coords[i] = i * DistStrip;
859  }
860 
861  //get and output moving body variables
862  int nStrVars = 3;
863  Array<OneD, Array<OneD, NekDouble> > Motion_x(nStrVars);
864  Array<OneD, Array<OneD, NekDouble> > Motion_y(nStrVars);
865 
866  for(int i = 0; i < nStrVars; i++)
867  {
868  Motion_x[i] = Array<OneD, NekDouble>(local_planes,0.0);
869  Motion_y[i] = Array<OneD, NekDouble>(local_planes,0.0);
870  }
871 
872  for(int plane = 0; plane < local_planes; plane++)
873  {
874  for (int var = 0; var < nStrVars; var++)
875  {
876  int xoffset = plane*nStrVars+var;
877  int yoffset = nStrVars*local_planes+xoffset;
878  Motion_x[var][plane] = MotionVars[xoffset];
879  Motion_y[var][plane] = MotionVars[yoffset];
880  }
881  }
882 
883  Array <OneD, NekDouble> CableMotions(6);
884 
885  for(int var = 0; var <nStrVars; var++)
886  {
887  CableMotions[var] = Motion_x[var][0];
888  CableMotions[3+var] = Motion_y[var][0];
889  }
890  // Send to root process.
891  if (colrank == 0)
892  {
893  m_outputStream[1].width(8);
894  m_outputStream[1] << setprecision(6) << time;
895  m_outputStream[1].width(25);
896  m_outputStream[1] << setprecision(6) << z_coords[0];
897  for(int var = 0; var < 2*nStrVars; var++)
898  {
899  m_outputStream[1].width(25);
900  m_outputStream[1] << setprecision(8) << CableMotions[var];
901  }
902  m_outputStream[1] << endl;
903 
904  for (int i = 1; i < nstrips; ++i)
905  {
906  vColComm->Recv(i, CableMotions);
907 
908  m_outputStream[1].width(8);
909  m_outputStream[1] << setprecision(6) << time;
910  m_outputStream[1].width(25);
911  m_outputStream[1] << setprecision(6) << z_coords[i];
912  for(int var = 0; var < 2*nStrVars; var++)
913  {
914  m_outputStream[1].width(25);
915  m_outputStream[1] << setprecision(8) << CableMotions[var];
916  }
917  m_outputStream[1] << endl;
918  }
919  }
920  else
921  {
922  for(int i = 1; i < nstrips; i++)
923  {
924  if(colrank == i)
925  {
926  vColComm->Send(0, CableMotions);
927  }
928  }
929  }
930  }
931 }
932 
933 
934 /**
935  *
936  */
937 void FilterMovingBody::v_Finalise(
938  const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields,
939  const NekDouble &time)
940 {
941  if (pFields[0]->GetComm()->GetRank() == 0)
942  {
943  m_outputStream[0].close();
944  m_outputStream[1].close();
945  }
946 }
947 
948 
949 /**
950  *
951  */
952 bool FilterMovingBody::v_IsTimeDependent()
953 {
954  return true;
955 }
956 }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:135
Nektar::FilterMovingBody::m_BoundaryString
std::string m_BoundaryString
Definition: FilterMovingBody.h:114
Nektar::FilterMovingBody::m_outputStream
Array< OneD, std::ofstream > m_outputStream
Definition: FilterMovingBody.h:117
Nektar::FilterMovingBody::v_IsTimeDependent
virtual bool v_IsTimeDependent()
Definition: FilterMovingBody.cpp:952
Expansion2D.h
Vmath::Vvtvp
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:428
Nektar::FilterMovingBody::m_outputFile_mot
std::string m_outputFile_mot
Definition: FilterMovingBody.h:119
Nektar::FilterMovingBody::FilterMovingBody
FilterMovingBody(const LibUtilities::SessionReaderSharedPtr &pSession, const std::map< std::string, std::string > &pParams)
Definition: FilterMovingBody.cpp:53
Nektar::FilterMovingBody::m_index
unsigned int m_index
Definition: FilterMovingBody.h:107
Nektar::FilterMovingBody::m_session
LibUtilities::SessionReaderSharedPtr m_session
Definition: FilterMovingBody.h:100
Nektar::LibUtilities::SessionReaderSharedPtr
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:50
Nektar::ParseUtils::GenerateSeqVector
static bool GenerateSeqVector(const char *const str, std::vector< unsigned int > &vec)
Definition: ParseUtils.hpp:78
Nektar::FilterMovingBody::v_Finalise
virtual void v_Finalise(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
Definition: FilterMovingBody.cpp:937
Nektar::FilterMovingBody::UpdateMotion
void UpdateMotion(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, Array< OneD, NekDouble > &MotionVars, const NekDouble &time)
Definition: FilterMovingBody.cpp:698
Nektar::StdRegions::StdExpansion::GetTotPoints
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
Definition: StdExpansion.h:141
Nektar::FilterMovingBody::create
static SolverUtils::FilterSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const std::map< std::string, std::string > &pParams)
Creates an instance of this class.
Definition: FilterMovingBody.h:58
Nektar::FilterMovingBody::className
static std::string className
Definition: FilterMovingBody.h:66
Nektar::FilterMovingBody::v_Initialise
virtual void v_Initialise(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
Definition: FilterMovingBody.cpp:124
Nektar::LibUtilities::CommSharedPtr
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:53
Nektar::FilterMovingBody::m_boundaryRegionsIdList
vector< unsigned int > m_boundaryRegionsIdList
ID's of boundary regions where we want the forces.
Definition: FilterMovingBody.h:103
Nektar::FilterMovingBody::m_isHomogeneous1D
bool m_isHomogeneous1D
Definition: FilterMovingBody.h:112
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:199
Nektar::FilterMovingBody::m_outputFrequency
unsigned int m_outputFrequency
Definition: FilterMovingBody.h:108
FilterMovingBody.h
Nektar::FilterMovingBody::~FilterMovingBody
~FilterMovingBody()
Definition: FilterMovingBody.cpp:115
NekMemoryManager.hpp
Nektar::SpatialDomains::BoundaryRegionCollection
std::map< int, BoundaryRegionShPtr > BoundaryRegionCollection
Definition: Conditions.h:204
npts
static std::string npts
Definition: InputFld.cpp:43
Nektar::FilterMovingBody::m_outputFile_fce
std::string m_outputFile_fce
Definition: FilterMovingBody.h:118
Expansion3D.h
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Vmath::Sadd
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha + x.
Definition: Vmath.cpp:301
Nektar::FilterMovingBody::m_boundaryRegionIsInList
vector< bool > m_boundaryRegionIsInList
Determines if a given Boundary Region is in m_boundaryRegionsIdList.
Definition: FilterMovingBody.h:106
Nektar::FilterMovingBody::UpdateForce
void UpdateForce(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, Array< OneD, NekDouble > &Aeroforces, const NekDouble &time)
Definition: FilterMovingBody.cpp:227
Nektar::StdRegions::find
InputIterator find(InputIterator first, InputIterator last, InputIterator startingpoint, const EqualityComparable &value)
Definition: StdRegions.hpp:312
Nektar::SolverUtils::GetFilterFactory
FilterFactory & GetFilterFactory()
Definition: Filter.cpp:42
Nektar::SpatialDomains::BoundaryConditions::GetBoundaryRegions
const BoundaryRegionCollection & GetBoundaryRegions(void) const
Definition: Conditions.h:225
Expansion1D.h
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
Nektar::StdRegions::StdExpansionSharedPtr
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
Definition: StdExpansion.h:1791
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1016
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:285
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:169
Generated on Tue Dec 8 2015 23:41:28 for Nektar++ by   doxygen 1.8.9.1