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