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// Permission is hereby granted, free of charge, to any person obtaining a
14// copy of this software and associated documentation files (the "Software"),
15// to deal in the Software without restriction, including without limitation
16// the rights to use, copy, modify, merge, publish, distribute, sublicense,
17// and/or sell copies of the Software, and to permit persons to whom the
18// Software is furnished to do so, subject to the following conditions:
19//
20// The above copyright notice and this permission notice shall be included
21// in all copies or substantial portions of the Software.
22//
23// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29// DEALINGS IN THE SOFTWARE.
30//
31// Description: Output moving body forces during time-stepping.
32//
33///////////////////////////////////////////////////////////////////////////////
34
35#include <IncNavierStokesSolver/Filters/FilterMovingBody.h>
36#include <LibUtilities/BasicUtils/ParseUtils.h>
37#include <LibUtilities/Memory/NekMemoryManager.hpp>
38#include <LocalRegions/Expansion1D.h>
39#include <LocalRegions/Expansion2D.h>
40#include <LocalRegions/Expansion3D.h>
41#include <iomanip>
42
43using namespace std;
44
45namespace Nektar
46{
47
48std::string FilterMovingBody::className =
49 SolverUtils::GetFilterFactory().RegisterCreatorFunction(
50 "MovingBody", FilterMovingBody::create, "Moving Body Filter");
51/**
52 *
53 */
54FilterMovingBody::FilterMovingBody(
55 const LibUtilities::SessionReaderSharedPtr &pSession,
56 const std::weak_ptr<SolverUtils::EquationSystem> &pEquation,
57 const ParamMap &pParams)
58 : Filter(pSession, pEquation)
59{
60 // OutputFile
61 auto it = pParams.find("OutputFile");
62 if (it == pParams.end())
63 {
64 m_outputFile_fce = pSession->GetSessionName();
65 m_outputFile_mot = pSession->GetSessionName();
66 }
67 else
68 {
69 ASSERTL0(it->second.length() > 0, "Missing parameter 'OutputFile'.");
70
71 m_outputFile_fce = it->second;
72 m_outputFile_mot = it->second;
73 }
74 if (!(m_outputFile_fce.length() >= 4 &&
75 m_outputFile_fce.substr(m_outputFile_fce.length() - 4) == ".fce"))
76 {
77 m_outputFile_fce += ".fce";
78 }
79 if (!(m_outputFile_mot.length() >= 4 &&
80 m_outputFile_mot.substr(m_outputFile_mot.length() - 4) == ".mot"))
81 {
82 m_outputFile_mot += ".mot";
83 }
84
85 // OutputFrequency
86 it = pParams.find("OutputFrequency");
87 if (it == pParams.end())
88 {
89 m_outputFrequency = 1;
90 }
91 else
92 {
93 LibUtilities::Equation equ(m_session->GetInterpreter(), it->second);
94 m_outputFrequency = round(equ.Evaluate());
95 }
96
97 pSession->MatchSolverInfo("Homogeneous", "1D", m_isHomogeneous1D, false);
98 ASSERTL0(m_isHomogeneous1D, "Moving Body implemented just for 3D "
99 "Homogeneous 1D discetisations.");
100
101 // Boundary (to calculate the forces)
102 it = pParams.find("Boundary");
103 ASSERTL0(it != pParams.end(), "Missing parameter 'Boundary'.");
104 ASSERTL0(it->second.length() > 0, "Missing parameter 'Boundary'.");
105 m_BoundaryString = it->second;
106}
107
108/**
109 *
110 */
111FilterMovingBody::~FilterMovingBody()
112{
113}
114
115/**
116 *
117 */
118void FilterMovingBody::v_Initialise(
119 const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields,
120 const NekDouble &time)
121{
122 boost::ignore_unused(time);
123
124 m_index_f = 0;
125 m_index_m = 0;
126 m_outputStream = Array<OneD, std::ofstream>(2);
127 // Parse the boundary regions into a list.
128 std::string::size_type FirstInd = m_BoundaryString.find_first_of('[') + 1;
129 std::string::size_type LastInd = m_BoundaryString.find_last_of(']') - 1;
130
131 ASSERTL0(FirstInd <= LastInd,
132 (std::string("Error reading boundary region definition:") +
133 m_BoundaryString)
134 .c_str());
135
136 std::string IndString =
137 m_BoundaryString.substr(FirstInd, LastInd - FirstInd + 1);
138
139 bool parseGood =
140 ParseUtils::GenerateSeqVector(IndString, m_boundaryRegionsIdList);
141
142 ASSERTL0(parseGood && !m_boundaryRegionsIdList.empty(),
143 (std::string("Unable to read boundary regions index "
144 "range for FilterAeroForces: ") +
145 IndString)
146 .c_str());
147
148 // determine what boundary regions need to be considered
149 size_t numBoundaryRegions = pFields[0]->GetBndConditions().size();
150
151 m_boundaryRegionIsInList.insert(m_boundaryRegionIsInList.end(),
152 numBoundaryRegions, 0);
153
154 SpatialDomains::BoundaryConditions bcs(m_session, pFields[0]->GetGraph());
155
156 const SpatialDomains::BoundaryRegionCollection &bregions =
157 bcs.GetBoundaryRegions();
158
159 size_t cnt = 0;
160 for (auto &it : bregions)
161 {
162 if (std::find(m_boundaryRegionsIdList.begin(),
163 m_boundaryRegionsIdList.end(),
164 it.first) != m_boundaryRegionsIdList.end())
165 {
166 m_boundaryRegionIsInList[cnt] = 1;
167 }
168 ++cnt;
169 }
170
171 LibUtilities::CommSharedPtr vComm = pFields[0]->GetComm();
172
173 if (vComm->GetRank() == 0)
174 {
175 // Open output stream for cable forces
176 m_outputStream[0].open(m_outputFile_fce.c_str());
177 m_outputStream[0] << "#";
178 m_outputStream[0].width(7);
179 m_outputStream[0] << "Time";
180 m_outputStream[0].width(15);
181 m_outputStream[0] << "z";
182 m_outputStream[0].width(15);
183 m_outputStream[0] << "Fx (press)";
184 m_outputStream[0].width(15);
185 m_outputStream[0] << "Fx (visc)";
186 m_outputStream[0].width(15);
187 m_outputStream[0] << "Fx (tot)";
188 m_outputStream[0].width(15);
189 m_outputStream[0] << "Fy (press)";
190 m_outputStream[0].width(15);
191 m_outputStream[0] << "Fy (visc)";
192 m_outputStream[0].width(15);
193 m_outputStream[0] << "Fy (tot)";
194 m_outputStream[0] << endl;
195
196 // Open output stream for cable motions
197 m_outputStream[1].open(m_outputFile_mot.c_str());
198 m_outputStream[1] << "#";
199 m_outputStream[1].width(7);
200 m_outputStream[1] << "Time";
201 m_outputStream[1].width(15);
202 m_outputStream[1] << "z";
203 m_outputStream[1].width(15);
204 m_outputStream[1] << "Disp_x";
205 m_outputStream[1].width(15);
206 m_outputStream[1] << "Vel_x";
207 m_outputStream[1].width(15);
208 m_outputStream[1] << "Acel_x";
209 m_outputStream[1].width(15);
210 m_outputStream[1] << "Disp_y";
211 m_outputStream[1].width(15);
212 m_outputStream[1] << "Vel_y";
213 m_outputStream[1].width(15);
214 m_outputStream[1] << "Acel_y";
215 m_outputStream[1] << endl;
216 }
217}
218
219/**
220 *
221 */
222void FilterMovingBody::UpdateForce(
223 const LibUtilities::SessionReaderSharedPtr &pSession,
224 const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields,
225 Array<OneD, NekDouble> &Aeroforces, const NekDouble &time)
226{
227 size_t n, cnt, elmtid, nq, offset, boundary;
228 size_t nt = pFields[0]->GetNpoints();
229 size_t dim = pFields.size() - 1;
230
231 LocalRegions::ExpansionSharedPtr elmt;
232 Array<OneD, int> BoundarytoElmtID;
233 Array<OneD, int> BoundarytoTraceID;
234 Array<OneD, MultiRegions::ExpListSharedPtr> BndExp;
235
236 Array<OneD, const NekDouble> P(nt);
237 Array<OneD, const NekDouble> U(nt);
238 Array<OneD, const NekDouble> V(nt);
239 Array<OneD, const NekDouble> W(nt);
240
241 Array<OneD, Array<OneD, NekDouble>> gradU(dim);
242 Array<OneD, Array<OneD, NekDouble>> gradV(dim);
243 Array<OneD, Array<OneD, NekDouble>> gradW(dim);
244
245 Array<OneD, Array<OneD, NekDouble>> fgradU(dim);
246 Array<OneD, Array<OneD, NekDouble>> fgradV(dim);
247 Array<OneD, Array<OneD, NekDouble>> fgradW(dim);
248
249 LibUtilities::CommSharedPtr vComm = pFields[0]->GetComm();
250 LibUtilities::CommSharedPtr vRowComm = vComm->GetRowComm();
251 LibUtilities::CommSharedPtr vColComm = vComm->GetColumnComm();
252 LibUtilities::CommSharedPtr vCColComm = vColComm->GetColumnComm();
253
254 // set up storage space for forces on all the planes (z-positions)
255 // on each processors some of them will remain empty as we may
256 // have just few planes per processor
257 size_t Num_z_pos = pFields[0]->GetHomogeneousBasis()->GetNumModes();
258 Array<OneD, NekDouble> Fx(Num_z_pos, 0.0);
259 Array<OneD, NekDouble> Fxp(Num_z_pos, 0.0);
260 Array<OneD, NekDouble> Fxv(Num_z_pos, 0.0);
261 Array<OneD, NekDouble> Fy(Num_z_pos, 0.0);
262 Array<OneD, NekDouble> Fyp(Num_z_pos, 0.0);
263 Array<OneD, NekDouble> Fyv(Num_z_pos, 0.0);
264
265 NekDouble rho = (pSession->DefinesParameter("rho"))
266 ? (pSession->GetParameter("rho"))
267 : 1;
268 NekDouble mu = rho * pSession->GetParameter("Kinvis");
269
270 for (size_t i = 0; i < pFields.size(); ++i)
271 {
272 pFields[i]->SetWaveSpace(false);
273 pFields[i]->BwdTrans(pFields[i]->GetCoeffs(), pFields[i]->UpdatePhys());
274 pFields[i]->SetPhysState(true);
275 }
276
277 // Get the number of local planes on the process and their IDs
278 // to properly locate the forces in the Fx, Fy etc. vectors.
279 Array<OneD, unsigned int> ZIDs;
280 ZIDs = pFields[0]->GetZIDs();
281 size_t local_planes = ZIDs.size();
282
283 // Homogeneous 1D case Compute forces on all WALL boundaries
284 // This only has to be done on the zero (mean) Fourier mode.
285 for (size_t plane = 0; plane < local_planes; plane++)
286 {
287 pFields[0]->GetPlane(plane)->GetBoundaryToElmtMap(BoundarytoElmtID,
288 BoundarytoTraceID);
289 BndExp = pFields[0]->GetPlane(plane)->GetBndCondExpansions();
290 StdRegions::StdExpansionSharedPtr bc;
291
292 // loop over the types of boundary conditions
293 for (cnt = n = 0; n < BndExp.size(); ++n)
294 {
295 if (m_boundaryRegionIsInList[n] == 1)
296 {
297 size_t nExp = BndExp[n]->GetExpSize();
298 for (size_t i = 0; i < nExp; ++i, cnt++)
299 {
300 // find element of this expansion.
301 elmtid = BoundarytoElmtID[cnt];
302 elmt = pFields[0]->GetPlane(plane)->GetExp(elmtid);
303 nq = elmt->GetTotPoints();
304 offset =
305 pFields[0]->GetPlane(plane)->GetPhys_Offset(elmtid);
306
307 // Initialise local arrays for the velocity
308 // gradients size of total number of quadrature
309 // points for each element (hence local).
310 for (size_t j = 0; j < dim; ++j)
311 {
312 gradU[j] = Array<OneD, NekDouble>(nq, 0.0);
313 gradV[j] = Array<OneD, NekDouble>(nq, 0.0);
314 gradW[j] = Array<OneD, NekDouble>(nq, 0.0);
315 }
316
317 // identify boundary of element
318 boundary = BoundarytoTraceID[cnt];
319
320 // Extract fields
321 U = pFields[0]->GetPlane(plane)->GetPhys() + offset;
322 V = pFields[1]->GetPlane(plane)->GetPhys() + offset;
323 P = pFields[3]->GetPlane(plane)->GetPhys() + offset;
324
325 // compute the gradients
326 elmt->PhysDeriv(U, gradU[0], gradU[1]);
327 elmt->PhysDeriv(V, gradV[0], gradV[1]);
328
329 // Get face 1D expansion from element expansion
330 bc = BndExp[n]->GetExp(i)->as<LocalRegions::Expansion1D>();
331
332 // number of points on the boundary
333 size_t nbc = bc->GetTotPoints();
334
335 // several vectors for computing the forces
336 Array<OneD, NekDouble> Pb(nbc, 0.0);
337
338 for (size_t j = 0; j < dim; ++j)
339 {
340 fgradU[j] = Array<OneD, NekDouble>(nbc, 0.0);
341 fgradV[j] = Array<OneD, NekDouble>(nbc, 0.0);
342 }
343
344 Array<OneD, NekDouble> drag_t(nbc, 0.0);
345 Array<OneD, NekDouble> lift_t(nbc, 0.0);
346 Array<OneD, NekDouble> drag_p(nbc, 0.0);
347 Array<OneD, NekDouble> lift_p(nbc, 0.0);
348 Array<OneD, NekDouble> temp(nbc, 0.0);
349 Array<OneD, NekDouble> temp2(nbc, 0.0);
350
351 // identify boundary of element .
352 boundary = BoundarytoTraceID[cnt];
353
354 // extraction of the pressure and wss on the
355 // boundary of the element
356 elmt->GetTracePhysVals(boundary, bc, P, Pb);
357
358 for (size_t j = 0; j < dim; ++j)
359 {
360 elmt->GetTracePhysVals(boundary, bc, gradU[j],
361 fgradU[j]);
362 elmt->GetTracePhysVals(boundary, bc, gradV[j],
363 fgradV[j]);
364 }
365
366 // normals of the element
367 const Array<OneD, Array<OneD, NekDouble>> &normals =
368 elmt->GetTraceNormal(boundary);
369
370 //
371 // Compute viscous tractive forces on wall from
372 //
373 // t_i = - T_ij * n_j (minus sign for force
374 // exerted BY fluid ON wall),
375 //
376 // where
377 //
378 // T_ij = viscous stress tensor (here in Cartesian
379 // coords)
380 // dU_i dU_j
381 // = RHO * KINVIS * ( ---- + ---- ) .
382 // dx_j dx_i
383
384 // a) DRAG TERMS
385 //-rho*kinvis*(2*du/dx*nx+(du/dy+dv/dx)*ny)
386
387 Vmath::Vadd(nbc, fgradU[1], 1, fgradV[0], 1, drag_t, 1);
388 Vmath::Vmul(nbc, drag_t, 1, normals[1], 1, drag_t, 1);
389
390 Vmath::Smul(nbc, 2.0, fgradU[0], 1, fgradU[0], 1);
391 Vmath::Vmul(nbc, fgradU[0], 1, normals[0], 1, temp2, 1);
392 Vmath::Smul(nbc, 0.5, fgradU[0], 1, fgradU[0], 1);
393
394 Vmath::Vadd(nbc, temp2, 1, drag_t, 1, drag_t, 1);
395 Vmath::Smul(nbc, -mu, drag_t, 1, drag_t, 1);
396
397 // zero temporary storage vector
398 Vmath::Zero(nbc, temp, 0.0);
399 Vmath::Zero(nbc, temp2, 0.0);
400
401 // b) LIFT TERMS
402 //-rho*kinvis*(2*dv/dy*ny+(du/dy+dv/dx)*nx)
403
404 Vmath::Vadd(nbc, fgradU[1], 1, fgradV[0], 1, lift_t, 1);
405 Vmath::Vmul(nbc, lift_t, 1, normals[0], 1, lift_t, 1);
406
407 Vmath::Smul(nbc, 2.0, fgradV[1], 1, fgradV[1], 1);
408 Vmath::Vmul(nbc, fgradV[1], 1, normals[1], 1, temp2, 1);
409 Vmath::Smul(nbc, -0.5, fgradV[1], 1, fgradV[1], 1);
410
411 Vmath::Vadd(nbc, temp2, 1, lift_t, 1, lift_t, 1);
412 Vmath::Smul(nbc, -mu, lift_t, 1, lift_t, 1);
413
414 // Compute normal tractive forces on all WALL
415 // boundaries
416
417 Vmath::Vvtvp(nbc, Pb, 1, normals[0], 1, drag_p, 1, drag_p,
418 1);
419 Vmath::Vvtvp(nbc, Pb, 1, normals[1], 1, lift_p, 1, lift_p,
420 1);
421
422 // integration over the boundary
423 Fxv[ZIDs[plane]] += bc->Integral(drag_t);
424 Fyv[ZIDs[plane]] += bc->Integral(lift_t);
425
426 Fxp[ZIDs[plane]] += bc->Integral(drag_p);
427 Fyp[ZIDs[plane]] += bc->Integral(lift_p);
428 }
429 }
430 else
431 {
432 cnt += BndExp[n]->GetExpSize();
433 }
434 }
435 }
436
437 for (size_t i = 0; i < pFields.size(); ++i)
438 {
439 pFields[i]->SetWaveSpace(true);
440 pFields[i]->BwdTrans(pFields[i]->GetCoeffs(), pFields[i]->UpdatePhys());
441 pFields[i]->SetPhysState(false);
442 }
443
444 // In case we are using an hybrid parallelisation we need
445 // to reduce the forces on the same plane coming from
446 // different mesh partitions.
447 // It is quite an expensive communication, therefore
448 // we check to make sure it is actually required.
449
450 if (vComm->GetRowComm()->GetSize() > 0)
451 {
452 // NOTE 1: We can eventually sum the viscous and pressure
453 // component before doing the communication, thus
454 // reducing by a factor 2 the communication.
455 // NOTE 2: We may want to set up in the Comm class an AllReduce
456 // routine wich can handle arrays more efficiently
457 for (size_t plane = 0; plane < local_planes; plane++)
458 {
459 vRowComm->AllReduce(Fxp[ZIDs[plane]], LibUtilities::ReduceSum);
460 vRowComm->AllReduce(Fxv[ZIDs[plane]], LibUtilities::ReduceSum);
461 vRowComm->AllReduce(Fyp[ZIDs[plane]], LibUtilities::ReduceSum);
462 vRowComm->AllReduce(Fyv[ZIDs[plane]], LibUtilities::ReduceSum);
463 }
464 }
465
466 // At this point on rank (0,n) of the Mesh partion communicator we have
467 // the total areo forces on the planes which are on the same column
468 // communicator. Since the planes are scattered on different processors
469 // some of the entries of the vector Fxp, Fxp etc. are still zero.
470 // Now we need to reduce the values on a single vector on rank (0,0) of the
471 // global communicator.
472 if (!pSession->DefinesSolverInfo("HomoStrip"))
473 {
474 if (vComm->GetRowComm()->GetRank() == 0)
475 {
476 for (size_t z = 0; z < Num_z_pos; z++)
477 {
478 vColComm->AllReduce(Fxp[z], LibUtilities::ReduceSum);
479 vColComm->AllReduce(Fxv[z], LibUtilities::ReduceSum);
480 vColComm->AllReduce(Fyp[z], LibUtilities::ReduceSum);
481 vColComm->AllReduce(Fyv[z], LibUtilities::ReduceSum);
482 }
483 }
484 }
485 else
486 {
487 if (vComm->GetRowComm()->GetRank() == 0)
488 {
489 for (size_t z = 0; z < Num_z_pos; z++)
490 {
491 vCColComm->AllReduce(Fxp[z], LibUtilities::ReduceSum);
492 vCColComm->AllReduce(Fxv[z], LibUtilities::ReduceSum);
493 vCColComm->AllReduce(Fyp[z], LibUtilities::ReduceSum);
494 vCColComm->AllReduce(Fyv[z], LibUtilities::ReduceSum);
495 }
496 }
497 }
498
499 if (!pSession->DefinesSolverInfo("HomoStrip"))
500 {
501 // set the forces imparted on the cable's wall
502 for (size_t plane = 0; plane < local_planes; plane++)
503 {
504 Aeroforces[plane] = Fxp[ZIDs[plane]] + Fxv[ZIDs[plane]];
505 Aeroforces[plane + local_planes] =
506 Fyp[ZIDs[plane]] + Fyv[ZIDs[plane]];
507 }
508
509 // Only output every m_outputFrequency.
510 if ((m_index_f++) % m_outputFrequency)
511 {
512 return;
513 }
514
515 // At thi point in rank (0,0) we have the full vectors
516 // containing Fxp,Fxv,Fyp and Fyv where different positions
517 // in the vectors correspond to different planes.
518 // Here we write it to file. We do it just on one porcess
519
520 Array<OneD, NekDouble> z_coords(Num_z_pos, 0.0);
521 Array<OneD, const NekDouble> pts =
522 pFields[0]->GetHomogeneousBasis()->GetZ();
523
524 NekDouble LZ;
525 pSession->LoadParameter("LZ", LZ);
526 Vmath::Smul(Num_z_pos, LZ / 2.0, pts, 1, z_coords, 1);
527 Vmath::Sadd(Num_z_pos, LZ / 2.0, z_coords, 1, z_coords, 1);
528 if (vComm->GetRank() == 0)
529 {
530
531 Vmath::Vadd(Num_z_pos, Fxp, 1, Fxv, 1, Fx, 1);
532 Vmath::Vadd(Num_z_pos, Fyp, 1, Fyv, 1, Fy, 1);
533
534 for (size_t i = 0; i < Num_z_pos; i++)
535 {
536 m_outputStream[0].width(8);
537 m_outputStream[0] << setprecision(6) << time;
538
539 m_outputStream[0].width(15);
540 m_outputStream[0] << setprecision(6) << z_coords[i];
541
542 m_outputStream[0].width(15);
543 m_outputStream[0] << setprecision(8) << Fxp[i];
544 m_outputStream[0].width(15);
545 m_outputStream[0] << setprecision(8) << Fxv[i];
546 m_outputStream[0].width(15);
547 m_outputStream[0] << setprecision(8) << Fx[i];
548
549 m_outputStream[0].width(15);
550 m_outputStream[0] << setprecision(8) << Fyp[i];
551 m_outputStream[0].width(15);
552 m_outputStream[0] << setprecision(8) << Fyv[i];
553 m_outputStream[0].width(15);
554 m_outputStream[0] << setprecision(8) << Fy[i];
555 m_outputStream[0] << endl;
556 }
557 }
558 }
559 else
560 {
561 // average the forces over local planes for each processor
562 Array<OneD, NekDouble> fces(6, 0.0);
563 for (size_t plane = 0; plane < local_planes; plane++)
564 {
565 fces[0] += Fxp[ZIDs[plane]] + Fxv[ZIDs[plane]];
566 fces[1] += Fyp[ZIDs[plane]] + Fyv[ZIDs[plane]];
567 fces[2] += Fxp[ZIDs[plane]];
568 fces[3] += Fyp[ZIDs[plane]];
569 fces[4] += Fxv[ZIDs[plane]];
570 fces[5] += Fyv[ZIDs[plane]];
571 }
572
573 fces[0] = fces[0] / local_planes;
574 fces[1] = fces[1] / local_planes;
575 fces[2] = fces[2] / local_planes;
576 fces[3] = fces[3] / local_planes;
577 fces[4] = fces[4] / local_planes;
578 fces[5] = fces[5] / local_planes;
579
580 // average the forces over communicators within each strip
581 vCColComm->AllReduce(fces[0], LibUtilities::ReduceSum);
582 vCColComm->AllReduce(fces[1], LibUtilities::ReduceSum);
583 vCColComm->AllReduce(fces[2], LibUtilities::ReduceSum);
584 vCColComm->AllReduce(fces[3], LibUtilities::ReduceSum);
585 vCColComm->AllReduce(fces[4], LibUtilities::ReduceSum);
586 vCColComm->AllReduce(fces[5], LibUtilities::ReduceSum);
587
588 size_t npts = vComm->GetColumnComm()->GetColumnComm()->GetSize();
589
590 fces[0] = fces[0] / npts;
591 fces[1] = fces[1] / npts;
592 fces[2] = fces[2] / npts;
593 fces[3] = fces[3] / npts;
594 fces[4] = fces[4] / npts;
595 fces[5] = fces[5] / npts;
596
597 for (size_t plane = 0; plane < local_planes; plane++)
598 {
599 Aeroforces[plane] = fces[0];
600 Aeroforces[plane + local_planes] = fces[1];
601 }
602
603 // Only output every m_outputFrequency.
604 if ((m_index_f) % m_outputFrequency)
605 {
606 return;
607 }
608
609 size_t colrank = vColComm->GetRank();
610 size_t nstrips;
611
612 NekDouble DistStrip;
613
614 pSession->LoadParameter("Strip_Z", nstrips);
615 pSession->LoadParameter("DistStrip", DistStrip);
616
617 Array<OneD, NekDouble> z_coords(nstrips);
618 for (size_t i = 0; i < nstrips; i++)
619 {
620 z_coords[i] = i * DistStrip;
621 }
622
623 if (colrank == 0)
624 {
625 m_outputStream[0].width(8);
626 m_outputStream[0] << setprecision(6) << time;
627
628 m_outputStream[0].width(15);
629 m_outputStream[0] << setprecision(6) << z_coords[0];
630
631 m_outputStream[0].width(15);
632 m_outputStream[0] << setprecision(8) << fces[2];
633 m_outputStream[0].width(15);
634 m_outputStream[0] << setprecision(8) << fces[4];
635 m_outputStream[0].width(15);
636 m_outputStream[0] << setprecision(8) << fces[0];
637
638 m_outputStream[0].width(15);
639 m_outputStream[0] << setprecision(8) << fces[3];
640 m_outputStream[0].width(15);
641 m_outputStream[0] << setprecision(8) << fces[5];
642 m_outputStream[0].width(15);
643 m_outputStream[0] << setprecision(8) << fces[1];
644 m_outputStream[0] << endl;
645
646 for (size_t i = 1; i < nstrips; i++)
647 {
648 vColComm->Recv(i, fces);
649
650 m_outputStream[0].width(8);
651 m_outputStream[0] << setprecision(6) << time;
652
653 m_outputStream[0].width(15);
654 m_outputStream[0] << setprecision(6) << z_coords[i];
655
656 m_outputStream[0].width(15);
657 m_outputStream[0] << setprecision(8) << fces[2];
658 m_outputStream[0].width(15);
659 m_outputStream[0] << setprecision(8) << fces[4];
660 m_outputStream[0].width(15);
661 m_outputStream[0] << setprecision(8) << fces[0];
662
663 m_outputStream[0].width(15);
664 m_outputStream[0] << setprecision(8) << fces[3];
665 m_outputStream[0].width(15);
666 m_outputStream[0] << setprecision(8) << fces[5];
667 m_outputStream[0].width(15);
668 m_outputStream[0] << setprecision(8) << fces[1];
669 m_outputStream[0] << endl;
670 }
671 }
672 else
673 {
674 for (size_t i = 1; i < nstrips; i++)
675 {
676 if (colrank == i)
677 {
678 vColComm->Send(0, fces);
679 }
680 }
681 }
682 }
683}
684
685/**
686 *
687 */
688void FilterMovingBody::UpdateMotion(
689 const LibUtilities::SessionReaderSharedPtr &pSession,
690 const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields,
691 Array<OneD, NekDouble> &MotionVars, const NekDouble &time)
692{
693 boost::ignore_unused(pFields);
694
695 // Only output every m_outputFrequency.
696 if ((m_index_m++) % m_outputFrequency)
697 {
698 return;
699 }
700
701 size_t npts;
702 // Length of the cable
703 NekDouble Length;
704
705 if (!pSession->DefinesSolverInfo("HomoStrip"))
706 {
707 pSession->LoadParameter("LZ", Length);
708 npts = m_session->GetParameter("HomModesZ");
709 }
710 else
711 {
712 pSession->LoadParameter("LC", Length);
713 npts = m_session->GetParameter("HomStructModesZ");
714 }
715
716 NekDouble z_coords;
717 for (size_t n = 0; n < npts; n++)
718 {
719 z_coords = Length / npts * n;
720 m_outputStream[1].width(8);
721 m_outputStream[1] << setprecision(6) << time;
722 m_outputStream[1].width(15);
723 m_outputStream[1] << setprecision(6) << z_coords;
724 m_outputStream[1].width(15);
725 m_outputStream[1] << setprecision(8) << MotionVars[n];
726 m_outputStream[1].width(15);
727 m_outputStream[1] << setprecision(8) << MotionVars[npts + n];
728 m_outputStream[1].width(15);
729 m_outputStream[1] << setprecision(8) << MotionVars[2 * npts + n];
730 m_outputStream[1].width(15);
731 m_outputStream[1] << setprecision(8) << MotionVars[3 * npts + n];
732 m_outputStream[1].width(15);
733 m_outputStream[1] << setprecision(8) << MotionVars[4 * npts + n];
734 m_outputStream[1].width(15);
735 m_outputStream[1] << setprecision(8) << MotionVars[5 * npts + n];
736 m_outputStream[1] << endl;
737 }
738}
739
740/**
741 *
742 */
743void FilterMovingBody::v_Finalise(
744 const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields,
745 const NekDouble &time)
746{
747 boost::ignore_unused(time);
748
749 if (pFields[0]->GetComm()->GetRank() == 0)
750 {
751 m_outputStream[0].close();
752 m_outputStream[1].close();
753 }
754}
755
756/**
757 *
758 */
759bool FilterMovingBody::v_IsTimeDependent()
760{
761 return true;
762}
763} // namespace Nektar
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
FilterMovingBody.h
Nektar::FilterMovingBody::m_outputFile_fce
std::string m_outputFile_fce
Definition: FilterMovingBody.h:120
Nektar::FilterMovingBody::m_outputFrequency
size_t m_outputFrequency
Definition: FilterMovingBody.h:110
Nektar::FilterMovingBody::m_outputFile_mot
std::string m_outputFile_mot
Definition: FilterMovingBody.h:121
Nektar::FilterMovingBody::v_Initialise
virtual void v_Initialise(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time) override
Definition: FilterMovingBody.cpp:118
Nektar::FilterMovingBody::className
static std::string className
Definition: FilterMovingBody.h:67
Nektar::FilterMovingBody::v_Finalise
virtual void v_Finalise(const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time) override
Definition: FilterMovingBody.cpp:743
Nektar::FilterMovingBody::m_BoundaryString
std::string m_BoundaryString
Definition: FilterMovingBody.h:116
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:688
Nektar::FilterMovingBody::v_IsTimeDependent
virtual bool v_IsTimeDependent() override
Definition: FilterMovingBody.cpp:759
Nektar::FilterMovingBody::m_boundaryRegionsIdList
std::vector< unsigned int > m_boundaryRegionsIdList
ID's of boundary regions where we want the forces.
Definition: FilterMovingBody.h:104
Nektar::FilterMovingBody::FilterMovingBody
FilterMovingBody(const LibUtilities::SessionReaderSharedPtr &pSession, const std::weak_ptr< SolverUtils::EquationSystem > &pEquation, const ParamMap &pParams)
Definition: FilterMovingBody.cpp:54
Nektar::FilterMovingBody::create
static SolverUtils::FilterSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const std::weak_ptr< SolverUtils::EquationSystem > &pEquation, const ParamMap &pParams)
Creates an instance of this class.
Definition: FilterMovingBody.h:56
Nektar::FilterMovingBody::m_index_f
size_t m_index_f
Definition: FilterMovingBody.h:108
Nektar::FilterMovingBody::m_outputStream
Array< OneD, std::ofstream > m_outputStream
Definition: FilterMovingBody.h:119
Nektar::FilterMovingBody::m_boundaryRegionIsInList
std::vector< bool > m_boundaryRegionIsInList
Determines if a given Boundary Region is in m_boundaryRegionsIdList.
Definition: FilterMovingBody.h:107
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:222
Nektar::FilterMovingBody::~FilterMovingBody
~FilterMovingBody()
Definition: FilterMovingBody.cpp:111
Nektar::FilterMovingBody::m_isHomogeneous1D
bool m_isHomogeneous1D
Definition: FilterMovingBody.h:114
Nektar::FilterMovingBody::m_index_m
size_t m_index_m
Definition: FilterMovingBody.h:109
Nektar::LibUtilities::NekFactory::RegisterCreatorFunction
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:198
Nektar::ParseUtils::GenerateSeqVector
static bool GenerateSeqVector(const std::string &str, std::vector< unsigned int > &out)
Takes a comma-separated compressed string and converts it to entries in a vector.
Definition: ParseUtils.cpp:105
Nektar::SolverUtils::Filter::m_session
LibUtilities::SessionReaderSharedPtr m_session
Definition: Filter.h:85
Nektar::SolverUtils::Filter::ParamMap
std::map< std::string, std::string > ParamMap
Definition: Filter.h:67
Nektar::SpatialDomains::BoundaryConditions::GetBoundaryRegions
const BoundaryRegionCollection & GetBoundaryRegions(void) const
Definition: Conditions.h:234
Nektar::StdRegions::StdExpansion::GetTotPoints
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
Definition: StdExpansion.h:140
Nektar::LibUtilities::SessionReaderSharedPtr
std::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: SessionReader.h:115
Nektar::LibUtilities::CommSharedPtr
std::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:57
Nektar::LibUtilities::P
@ P
Monomial polynomials .
Definition: BasisType.h:64
Nektar::LocalRegions::ExpansionSharedPtr
std::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
Nektar::SolverUtils::GetFilterFactory
FilterFactory & GetFilterFactory()
Definition: Filter.cpp:41
Nektar::SpatialDomains::BoundaryRegionCollection
std::map< int, BoundaryRegionShPtr > BoundaryRegionCollection
Definition: Conditions.h:210
Nektar::StdRegions::find
InputIterator find(InputIterator first, InputIterator last, InputIterator startingpoint, const EqualityComparable &value)
Definition: StdRegions.hpp:453
Nektar::StdRegions::StdExpansionSharedPtr
std::shared_ptr< StdExpansion > StdExpansionSharedPtr
Definition: StdExpansion.h:1775
Nektar::UnitTests::z
std::vector< double > z(NPUPPER)
Nektar
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:2
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
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:207
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:569
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:354
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.cpp:245
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:487
Vmath::Sadd
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add scalar y = alpha + x.
Definition: Vmath.cpp:379
Generated on Wed Jan 31 2024 08:26:50 for Nektar++ by doxygen 1.9.4