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