Nektar++
CouplingCwipi.cpp
Go to the documentation of this file.
1////////////////////////////////////////////////////////////////////////////////
2//
3// File: CouplingCwipi.cpp
4//
5// For more information, please see: http://www.nektar.info/
6//
7// The MIT License
8//
9// Copyright (c) 2017 Kilian Lackhove
10//
11// Permission is hereby granted, free of charge, to any person obtaining a
12// copy of this software and associated documentation files (the "Software"),
13// to deal in the Software without restriction, including without limitation
14// the rights to use, copy, modify, merge, publish, distribute, sublicense,
15// and/or sell copies of the Software, and to permit persons to whom the
16// Software is furnished to do so, subject to the following conditions:
17//
18// The above copyright notice and this permission notice shall be included
19// in all copies or substantial portions of the Software.
20//
21// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
22// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
23// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
24// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
25// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
26// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
27// DEALINGS IN THE SOFTWARE.
28//
29// Description: CWIPI Exchange class
30//
31////////////////////////////////////////////////////////////////////////////////
32
33#include "CouplingCwipi.h"
34
43
45
46#include <boost/format.hpp>
47
48#define OUTPUT_FREQ 0
49
50namespace Nektar::SolverUtils
51{
52
53std::string CouplingCwipi::className =
54 GetCouplingFactory().RegisterCreatorFunction("Cwipi", CouplingCwipi::create,
55 "Cwipi Coupling");
56
58 const int entities_dim, const int n_local_vertex,
59 [[maybe_unused]] const int n_local_element,
60 [[maybe_unused]] const int n_local_polyhedra, const int n_distant_point,
61 [[maybe_unused]] const double local_coordinates[],
62 [[maybe_unused]] const int local_connectivity_index[],
63 [[maybe_unused]] const int local_connectivity[],
64 [[maybe_unused]] const int local_polyhedra_face_index[],
65 [[maybe_unused]] const int local_polyhedra_cell_to_face_connectivity[],
66 [[maybe_unused]] const int local_polyhedra_face_connectivity_index[],
67 [[maybe_unused]] const int local_polyhedra_face_connectivity[],
68 const double distant_points_coordinates[],
69 [[maybe_unused]] const int distant_points_location[],
70 [[maybe_unused]] const float distant_points_distance[],
71 [[maybe_unused]] const int distant_points_barycentric_coordinates_index[],
72 [[maybe_unused]] const double distant_points_barycentric_coordinates[],
73 const int stride, [[maybe_unused]] const cwipi_solver_type_t solver_type,
74 [[maybe_unused]] const void *local_field, void *distant_field)
75{
76 Array<OneD, Array<OneD, NekDouble>> interpField(stride);
77
78 Array<OneD, Array<OneD, NekDouble>> distCoords(n_distant_point);
79 for (int i = 0; i < n_distant_point; ++i)
80 {
81 distCoords[i] = Array<OneD, NekDouble>(3);
82 for (int j = 0; j < 3; ++j)
83 {
84 distCoords[i][j] = distant_points_coordinates[3 * i + j];
85 }
86 }
87
88 std::stringstream sst;
89 sst << entities_dim << "," << n_local_vertex << "," << stride;
90 SendCallbackMap[sst.str()](interpField, distCoords);
91
92 ASSERTL0(interpField.size() == stride, "size mismatch");
93 ASSERTL0(interpField[0].size() == n_distant_point, "size mismatch");
94
95 for (int i = 0; i < n_distant_point; i++)
96 {
97 for (int j = 0; j < stride; ++j)
98 {
99 ((double *)distant_field)[i * stride + j] = interpField[j][i];
100 }
101 }
102}
103
105 : Coupling(field), m_sendHandle(-1), m_recvHandle(-1), m_lastSend(-1E6),
106 m_lastReceive(-1E6), m_points(nullptr), m_coords(nullptr),
107 m_connecIdx(nullptr), m_connec(nullptr), m_rValsInterl(nullptr),
108 m_sValsInterl(nullptr)
109{
110 // defaults
111 m_config["GEOMTOL"] = "0.1";
112 m_config["LOCALNAME"] = "nektar";
113 m_config["REMOTENAME"] = "precise";
114 m_config["OVERSAMPLE"] = "0";
115 m_config["FILTERWIDTH"] = "-1";
116 m_config["DUMPRAW"] = "0";
117 m_config["SENDMETHOD"] = "EVALUATE";
118 m_config["NOTLOCMETHOD"] = "KEEP";
119}
120
122{
124
125 ReadConfig(m_evalField->GetSession());
126
127 cwipi_add_local_int_control_parameter("nSendVars", m_nSendVars);
128 cwipi_add_local_int_control_parameter("nRecvVars", m_nRecvVars);
129 cwipi_add_local_string_control_parameter(
130 "recvFieldNames", m_config["RECEIVEVARIABLES"].c_str());
131 cwipi_add_local_string_control_parameter("sendFieldNames",
132 m_config["SENDVARIABLES"].c_str());
133
134 // MPI_TAG_UB is guaranteed to be at least 32767, so we make
135 // sure m_recvTag < 32767. Only caveat: m_recvTag is not guaranteed to be
136 // unique.
137 m_recvTag =
138 std::hash<std::string>()(m_couplingName + m_config["REMOTENAME"] +
139 m_config["LOCALNAME"]) %
140 32767;
141
142 cwipi_add_local_int_control_parameter("receiveTag", m_recvTag);
143
144 m_spacedim = m_evalField->GetGraph()->GetSpaceDimension();
145
146 if (m_filtWidth > 0)
147 {
148 m_filtWidth = 2 * M_PI / m_filtWidth;
150 }
151
152 // Init Coupling
153 cwipi_solver_type_t solver_type = CWIPI_SOLVER_CELL_VERTEX;
154 NekDouble geom_tol = std::stod(m_config["GEOMTOL"]);
155 cwipi_create_coupling(
156 m_couplingName.c_str(), CWIPI_COUPLING_PARALLEL_WITH_PARTITIONING,
157 m_config["REMOTENAME"].c_str(), m_spacedim, geom_tol, CWIPI_STATIC_MESH,
158 solver_type, OUTPUT_FREQ, "Ensight Gold", "text");
159 cwipi_synchronize_control_parameter(m_config["REMOTENAME"].c_str());
160
161 if (m_evalField->GetComm()->GetRank() == 0 &&
162 m_evalField->GetSession()->DefinesCmdLineArgument("verbose"))
163 {
164 cwipi_dump_application_properties();
165 }
166
167 m_sendTag = cwipi_get_distant_int_control_parameter(
168 m_config["REMOTENAME"].c_str(), "receiveTag");
169
170 if (cwipi_has_int_parameter(m_config["REMOTENAME"].c_str(), "nRecvVars"))
171 {
172 int remoteNRecvVars = cwipi_get_distant_int_control_parameter(
173 m_config["REMOTENAME"].c_str(), "nRecvVars");
174 ASSERTL0(remoteNRecvVars == m_nSendVars,
175 "Number of local send vars different to remote received vars");
176 }
177
178 if (cwipi_has_int_parameter(m_config["REMOTENAME"].c_str(), "nSendVars"))
179 {
180 int remoteNSendVars = cwipi_get_distant_int_control_parameter(
181 m_config["REMOTENAME"].c_str(), "nSendVars");
182 ASSERTL0(remoteNSendVars == m_nRecvVars,
183 "Number of local receive vars different to remote sent vars");
184 }
185
186 AnnounceMesh();
187
188 if (m_nRecvVars > 0 && m_recvSteps > 0)
189 {
190 SetupReceive();
191 }
192
193 if (m_evalField->GetComm()->GetRank() == 0 &&
194 m_evalField->GetSession()->DefinesCmdLineArgument("verbose"))
195 {
196 std::cout << "locating..." << std::endl;
197 }
198 cwipi_locate(m_couplingName.c_str());
199
200 if (m_nSendVars > 0 && m_sendSteps > 0)
201 {
202 SetupSend();
203 }
204
205 if (m_nRecvVars > 0 && m_recvSteps > 0)
206 {
207 ReceiveStart();
208 }
209}
210
212{
213 free(m_coords);
214 free(m_points);
215 free(m_connec);
216 free(m_connecIdx);
217 free(m_rValsInterl);
218 free(m_sValsInterl);
219}
220
222{
223 ASSERTL0(session->DefinesElement("Nektar/Coupling"),
224 "No Coupling config found");
225
226 m_config["LOCALNAME"] = session->GetCmdLineArgument<std::string>("cwipi");
227
228 TiXmlElement *vCoupling = session->GetElement("Nektar/Coupling");
229 ASSERTL0(vCoupling, "Invalid Coupling config");
230
231 m_filtWidth = std::stod(m_config["FILTERWIDTH"]);
232}
233
235{
236 int oversamp = std::stoi(m_config["OVERSAMPLE"]);
237
240 recvGraph->SetExpansionInfoToPointOrder(
241 oversamp + m_evalField->GetExp(0)->GetNumPoints(0));
242
243 // TODO: DeclareCoeffPhysArrays
245 m_evalField->GetSession(), recvGraph, "DefaultVar");
246
249 for (int i = 0; i < m_nRecvVars; ++i)
250 {
251 m_oldFields[i] = Array<OneD, NekDouble>(m_evalField->GetTotPoints());
252 m_newFields[i] = Array<OneD, NekDouble>(m_evalField->GetTotPoints());
253 }
254
255 // define the quadrature points at which we want to receive data
256 m_nPoints = m_recvField->GetTotPoints();
261 m_recvField->GetCoords(coords[0], coords[1], coords[2]);
262
263 m_points = (double *)malloc(sizeof(double) * 3 * m_nPoints);
264 ASSERTL1(m_points != nullptr, "malloc failed for m_points");
265
266 for (int i = 0; i < m_nPoints; ++i)
267 {
268 m_points[3 * i + 0] = double(coords[0][i]);
269
270 if (m_spacedim > 1)
271 {
272 m_points[3 * i + 1] = double(coords[1][i]);
273 }
274 else
275 {
276 m_points[3 * i + 1] = 0.0;
277 }
278
279 if (m_spacedim > 2)
280 {
281 m_points[3 * i + 2] = double(coords[2][i]);
282 }
283 else
284 {
285 m_points[3 * i + 2] = 0.0;
286 }
287 }
288
289 cwipi_set_points_to_locate(m_couplingName.c_str(), m_nPoints, m_points);
290
291 m_rValsInterl = (double *)malloc(sizeof(double) * m_nPoints * m_nRecvVars);
292 ASSERTL1(m_rValsInterl != nullptr, "malloc failed for m_rValsInterl");
293}
294
296{
297 // this array is never used because of our send callback method
298 m_sValsInterl = (double *)malloc(
299 sizeof(double) * m_evalField->GetGraph()->GetNvertices() * m_nSendVars);
300 ASSERTL1(m_sValsInterl != nullptr, "malloc failed for m_sValsInterl");
301 for (int i = 0; i < m_evalField->GetGraph()->GetNvertices() * m_nSendVars;
302 ++i)
303 {
304 m_sValsInterl[i] = i;
305 }
306
307 // register this coupling as sender
308 std::stringstream sst;
309 sst << m_spacedim << "," << m_evalField->GetGraph()->GetNvertices() << ","
310 << m_nSendVars;
311 SendCallbackMap[sst.str()] =
312 std::bind(&CouplingCwipi::SendCallback, this, std::placeholders::_1,
313 std::placeholders::_2);
314 cwipi_set_interpolation_function(m_couplingName.c_str(),
316}
317
319 Array<OneD, Array<OneD, NekDouble>> interpField,
320 Array<OneD, Array<OneD, NekDouble>> distCoords)
321{
322 int nOutPts = distCoords.size();
323
325 for (int i = 0; i < nOutPts; ++i)
326 {
327 // Obtain Element and LocalCoordinate to interpolate
328 int elmtid = -1;
330 while (elmtid < 0 && tol <= 1E3 * NekConstants::kNekZeroTol)
331 {
332 elmtid = m_evalField->GetExpIndex(distCoords[i], Lcoords, tol);
333 tol *= 2;
334 }
335 if (tol > 2 * NekConstants::kNekZeroTol)
336 {
337 for (int j = 0; j < m_spacedim; ++j)
338 {
339 if (Lcoords[j] < -1 - 0.75 * NekConstants::kNekZeroTol)
340 {
341 Lcoords[j] = -1;
342 }
343 if (Lcoords[j] > 1 + 0.75 * NekConstants::kNekZeroTol)
344 {
345 Lcoords[j] = 1;
346 }
347 }
348 }
349
350 ASSERTL0(elmtid >= 0,
351 "no element found for (" +
352 boost::lexical_cast<std::string>(distCoords[i][0]) + ", " +
353 boost::lexical_cast<std::string>(distCoords[i][1]) + ", " +
354 boost::lexical_cast<std::string>(distCoords[i][2]) + ")");
355
356 int offset = m_evalField->GetPhys_Offset(elmtid);
357
358 for (int f = 0; f < m_nSendVars; ++f)
359 {
360 NekDouble value = m_evalField->GetExp(elmtid)->StdPhysEvaluate(
361 Lcoords, m_sendField[f] + offset);
362
363 ASSERTL0(!(boost::math::isnan)(value), "new value is not a number");
364 interpField[f][i] = value;
365 }
366 }
367}
368
370{
371 const double *distCoords =
372 cwipi_get_distant_coordinates(m_couplingName.c_str());
373 int nPts = cwipi_get_n_distant_points(m_couplingName.c_str());
374
376 for (int i = 0; i < 3; ++i)
377 {
378 local[i] = Array<OneD, NekDouble>(m_evalField->GetTotPoints(), 0.0);
379 }
380 m_evalField->GetCoords(local[0], local[1], local[2]);
383
385 for (int i = 0; i < 3; ++i)
386 {
387 dist[i] = Array<OneD, NekDouble>(nPts);
388 for (int j = 0; j < nPts; ++j)
389 {
390 dist[i][j] = distCoords[3 * j + i];
391 }
392 }
395
397 if (boost::to_upper_copy(m_config["SENDMETHOD"]) == "SHEPARD")
398 {
399 method = LibUtilities::eShepard;
400 }
402 MultiRegions::ExpListSharedPtr>>>::AllocateSharedPtr(method);
403 m_sendInterpolator->CalcWeights(locatPts, distPts);
404 m_sendInterpolator->PrintStatistics();
405}
406
408{
410
411 // get Elements
419 if (m_spacedim == 1)
420 {
421 seggeom = graph->GetAllSegGeoms();
422 }
423 else if (m_spacedim == 2)
424 {
425 trigeom = graph->GetAllTriGeoms();
426 quadgeom = graph->GetAllQuadGeoms();
427 }
428 else if (m_spacedim == 3)
429 {
430 tetgeom = graph->GetAllTetGeoms();
431 pyrgeom = graph->GetAllPyrGeoms();
432 prismgeom = graph->GetAllPrismGeoms();
433 hexgeom = graph->GetAllHexGeoms();
434 };
435
436 int nVerts = graph->GetNvertices();
437 int nElts = seggeom.size() + trigeom.size() + quadgeom.size() +
438 tetgeom.size() + pyrgeom.size() + prismgeom.size() +
439 hexgeom.size();
440
441 // allocate CWIPI arrays
442 m_coords = (double *)malloc(sizeof(double) * 3 * nVerts);
443 ASSERTL1(m_coords != nullptr, "malloc failed for m_coords");
444 int tmp = 2 * seggeom.size() + 3 * trigeom.size() + 4 * quadgeom.size() +
445 4 * tetgeom.size() + 5 * pyrgeom.size() + 6 * prismgeom.size() +
446 8 * hexgeom.size();
447 m_connec = (int *)malloc(sizeof(int) * tmp);
448 ASSERTL1(m_connec != nullptr, "malloc failed for m_connec");
449 m_connecIdx = (int *)malloc(sizeof(int) * (nElts + 1));
450 ASSERTL1(m_connecIdx != nullptr, "malloc failed for m_connecIdx");
451
452 m_connecIdx[0] = 0;
453 int coordsPos = 0;
454 int connecPos = 0;
455 int conidxPos = 0;
456
457 AddElementsToMesh(seggeom, coordsPos, connecPos, conidxPos);
458 AddElementsToMesh(trigeom, coordsPos, connecPos, conidxPos);
459 AddElementsToMesh(quadgeom, coordsPos, connecPos, conidxPos);
460 AddElementsToMesh(tetgeom, coordsPos, connecPos, conidxPos);
461 AddElementsToMesh(pyrgeom, coordsPos, connecPos, conidxPos);
462 AddElementsToMesh(prismgeom, coordsPos, connecPos, conidxPos);
463 AddElementsToMesh(hexgeom, coordsPos, connecPos, conidxPos);
464
465 // output the mesh in tecplot format. If this works, CWIPI will be able
466 // to process it, too
467 /*
468 std::cout << "VARIABLES = \"X\", \"Y\", \"Z\", \"U\"" << std::endl;
469 std::cout << "ZONE NODES=" << nVerts << ",ELEMENTS=" << nElts;
470 std::cout << ",DATAPACKING=POINT, ZONETYPE=FETETRAHEDRON" << std::endl;
471 for (int i = 0; i < nVerts; ++i)
472 {
473 std::cout << m_coords[3*i + 0] << " " << m_coords[3*i + 1] << " ";
474 std::cout << m_coords[3*i + 2] << " " << 1.0 << std::endl;
475 }
476 for (int i = 0; i < nElts; ++i)
477 {
478 std::cout << m_connec[i*4 + 0] << " " << m_connec[i*4 + 1] << " ";
479 std::cout << m_connec[i*4 + 2] << " " << m_connec[i*4 + 3] <<
480 std::endl;
481 }
482 */
483
484 cwipi_define_mesh(m_couplingName.c_str(), nVerts, nElts, m_coords,
486}
487
489{
490 cwipi_delete_coupling(m_couplingName.c_str());
491}
492
493template <typename T>
494void CouplingCwipi::AddElementsToMesh(T geom, int &coordsPos, int &connecPos,
495 int &conidxPos)
496{
497 // helper variables
500 int vertID;
501
502 int kNverts = T::mapped_type::element_type::kNverts;
503
504 // iterate over all elements
505 for (auto it = geom.begin(); it != geom.end(); it++)
506 {
507 // iterate over the elements vertices
508 for (int j = 0; j < kNverts; ++j)
509 {
510 vert = it->second->GetVertex(j);
511 vertID = vert->GetVid();
512
513 // check if we already stored the vertex
514 if (m_vertMap.count(vertID) == 0)
515 {
516 // store the vertex
517 vert->GetCoords(x[0], x[1], x[2]);
518 m_coords[3 * coordsPos + 0] = double(x[0]);
519 m_coords[3 * coordsPos + 1] = double(x[1]);
520 m_coords[3 * coordsPos + 2] = double(x[2]);
521
522 // store the vertex position in the m_coords array
523 m_vertMap[vertID] = coordsPos;
524 coordsPos++;
525 }
526
527 m_connec[connecPos] = m_vertMap[vertID] + 1;
528 connecPos++;
529 }
530
531 m_connecIdx[conidxPos + 1] = m_connecIdx[conidxPos] + kNverts;
532 conidxPos++;
533 }
534}
535
537 Array<OneD, Array<OneD, NekDouble>> &interpField,
538 Array<OneD, Array<OneD, NekDouble>> &distCoords)
539{
540 ASSERTL0(interpField.size() == m_nSendVars, "size mismatch");
541
542 for (int i = 0; i < m_nSendVars; ++i)
543 {
544 interpField[i] = Array<OneD, NekDouble>(distCoords.size());
545 }
546
547 if (boost::to_upper_copy(m_config["SENDMETHOD"]) == "NEARESTNEIGHBOUR" ||
548 boost::to_upper_copy(m_config["SENDMETHOD"]) == "SHEPARD")
549 {
551 {
553 }
554
557 0, m_sendField);
560 0, interpField);
561 m_sendInterpolator->Interpolate(ptsIn, ptsOut);
562 }
563 else
564 {
565 EvaluateFields(interpField, distCoords);
566 }
567}
568
570 const int step, const NekDouble time,
571 const Array<OneD, const Array<OneD, NekDouble>> &field,
572 std::vector<std::string> &varNames)
573{
574 if (m_nSendVars < 1 || m_sendSteps < 1)
575 {
576 return;
577 }
578
579 if (step >= m_lastSend + m_sendSteps)
580 {
581 SendComplete();
582
583 m_lastSend = step;
584
585 std::vector<int> sendVarsToVars =
588 for (int i = 0; i < sendVarsToVars.size(); ++i)
589 {
590 m_sendField[i] = field[sendVarsToVars[i]];
591 }
592
593 if (m_evalField->GetComm()->GetRank() == 0 &&
594 m_evalField->GetSession()->DefinesCmdLineArgument("verbose"))
595 {
596 std::cout << "sending fields at i = " << step << ", t = " << time
597 << std::endl;
598 }
599
600 LibUtilities::Timer timer1;
601 timer1.Start();
602
603 char sendFN[10];
604 strcpy(sendFN, "dummyName");
605
606 cwipi_issend(m_couplingName.c_str(), "ex1", m_sendTag, m_nSendVars,
607 step, time, sendFN, m_sValsInterl, &m_sendHandle);
608 timer1.Stop();
609
610 if (m_evalField->GetComm()->GetRank() == 0 &&
611 m_evalField->GetSession()->DefinesCmdLineArgument("verbose"))
612 {
613 std::cout << "Send total time: " << timer1.TimePerTest(1)
614 << std::endl;
615 }
616 }
617}
618
620{
621 if (m_sendHandle < 0)
622 {
623 return;
624 }
625
626 LibUtilities::Timer timer1;
627 timer1.Start();
628 cwipi_wait_issend(m_couplingName.c_str(), m_sendHandle);
629 timer1.Stop();
630
631 // set to -1 so we dont try finishing a send before a new one was started
632 m_sendHandle = -1;
633
634 if (m_evalField->GetComm()->GetRank() == 0 &&
635 m_evalField->GetSession()->DefinesCmdLineArgument("verbose"))
636 {
637 std::cout << "Send waiting time: " << timer1.TimePerTest(1)
638 << std::endl;
639 }
640}
641
643{
644 if (m_recvHandle >= 0)
645 {
646 return;
647 }
648
649 LibUtilities::Timer timer1;
650 timer1.Start();
651 // workaround a bug in cwipi: receiving_field_name should be const char* but
652 // is char*
653 char recFN[10];
654 strcpy(recFN, "dummyName");
655
656 cwipi_irecv(m_couplingName.c_str(), "ex1", m_recvTag, m_nRecvVars, 0, 0.0,
657 recFN, m_rValsInterl, &m_recvHandle);
658 timer1.Stop();
659
660 if (m_evalField->GetComm()->GetRank() == 0 &&
661 m_evalField->GetSession()->DefinesCmdLineArgument("verbose"))
662 {
663 std::cout << "Receive start time: " << timer1.TimePerTest(1)
664 << std::endl;
665 }
666}
667
668void CouplingCwipi::v_Receive(const int step, const NekDouble time,
670 std::vector<std::string> &varNames)
671{
672 if (m_nRecvVars < 1 || m_recvSteps < 1)
673 {
674 return;
675 }
676
678 std::vector<int> recvVarsToVars =
680 ASSERTL1(m_nRecvVars == recvVarsToVars.size(), "field size mismatch");
681 for (int i = 0; i < recvVarsToVars.size(); ++i)
682 {
683 recvFields[i] = field[recvVarsToVars[i]];
684 }
685
686 int nq = m_evalField->GetTotPoints();
687
688 // make sure we have sensible data in old/new recvFields the first time this
689 // method is called
690 if (m_lastReceive < 0)
691 {
692 for (int i = 0; i < m_nRecvVars; ++i)
693 {
694 Vmath::Vcopy(nq, recvFields[i], 1, m_oldFields[i], 1);
695 Vmath::Vcopy(nq, recvFields[i], 1, m_newFields[i], 1);
696 }
697 }
698
699 if (step >= m_lastReceive + m_recvSteps)
700 {
701 for (int i = 0; i < m_nRecvVars; ++i)
702 {
703 Vmath::Vcopy(nq, m_newFields[i], 1, m_oldFields[i], 1);
704 }
705
706 ReceiveCwipi(step, time, m_newFields);
707 }
708
709 NekDouble fact =
711 for (int i = 0; i < m_nRecvVars; ++i)
712 {
713 Vmath::Svtsvtp(nq, fact, m_newFields[i], 1, (1 - fact), m_oldFields[i],
714 1, recvFields[i], 1);
715 }
716}
717
718void CouplingCwipi::ReceiveCwipi(const int step, const NekDouble time,
720{
721 ASSERTL1(m_nRecvVars == field.size(), "field size mismatch");
722
723 if (m_nRecvVars < 1 || m_recvSteps < 1)
724 {
725 return;
726 }
727
728 if (step >= m_lastReceive + m_recvSteps)
729 {
730 m_lastReceive = step;
731
732 if (m_evalField->GetComm()->GetRank() == 0 &&
733 m_evalField->GetSession()->DefinesCmdLineArgument("verbose"))
734 {
735 std::cout << "waiting for receive at i = " << step
736 << ", t = " << time << std::endl;
737 }
738
739 LibUtilities::Timer timer1, timer2, timer3;
740 timer1.Start();
741
742 Array<OneD, NekDouble> tmpC(m_recvField->GetNcoeffs());
744 for (int i = 0; i < m_nRecvVars; ++i)
745 {
746 rVals[i] = Array<OneD, NekDouble>(m_recvField->GetTotPoints());
747 }
748
749 timer2.Start();
750 cwipi_wait_irecv(m_couplingName.c_str(), m_recvHandle);
751 timer2.Stop();
752
753 // set to -1 so we know we can start receiving again
754 m_recvHandle = -1;
755
756 int nNotLoc = cwipi_get_n_not_located_points(m_couplingName.c_str());
757 Array<OneD, int> notLoc;
758 if (nNotLoc != 0)
759 {
760 std::cout << "WARNING: relocating " << nNotLoc << " of "
761 << m_nPoints << " points" << std::endl;
762
763 const int *tmp =
764 cwipi_get_not_located_points(m_couplingName.c_str());
765 notLoc = Array<OneD, int>(nNotLoc);
766 for (int i = 0; i < nNotLoc; ++i)
767 {
768 notLoc[i] = tmp[i] - 1;
769 }
770
771 if (boost::to_upper_copy(m_config["NOTLOCMETHOD"]) == "KEEP")
772 {
773 // interpolate from m_evalField to m_recvField
774 for (int i = 0; i < m_nRecvVars; ++i)
775 {
776 m_evalField->FwdTrans(field[i], tmpC);
777 m_recvField->BwdTrans(tmpC, rVals[i]);
778 }
779 }
780 }
781
782 for (int i = 0, locPos = 0, intPos = 0; i < m_nPoints; ++i)
783 {
784 if (locPos < nNotLoc && notLoc[locPos] == i)
785 {
786 // keep the original value of field[j][i]
787 locPos++;
788 }
789 else
790 {
791 for (int j = 0; j < m_nRecvVars; ++j)
792 {
793 rVals[j][i] = m_rValsInterl[intPos * m_nRecvVars + j];
794 }
795 intPos++;
796 }
797 }
798
799 if (boost::to_upper_copy(m_config["NOTLOCMETHOD"]) == "EXTRAPOLATE")
800 {
801 int doExtrapolate = 0;
802 if (nNotLoc != 0)
803 {
804 doExtrapolate = 1;
805 }
806 m_evalField->GetSession()->GetComm()->AllReduce(
807 doExtrapolate, LibUtilities::ReduceMax);
808 if (doExtrapolate > 0)
809 {
810 ExtrapolateFields(rVals, notLoc);
811 }
812 }
813
814 if (m_config["DUMPRAW"] != "0")
815 {
816 DumpRawFields(time, rVals);
817 }
818
819 if (m_filtWidth > 0)
820 {
821 for (int i = 0; i < m_nRecvVars; ++i)
822 {
823 timer3.Start();
824
826
828 for (int j = 0; j < m_spacedim; ++j)
829 {
830 Velocity[j] = Array<OneD, NekDouble>(m_nPoints, 0.0);
831 }
832
833 Vmath::Smul(m_nPoints, -m_filtWidth, rVals[i], 1, forcing, 1);
834
836 StdRegions::VarCoeffMap varcoeffs;
837 MultiRegions::VarFactorsMap varfactors =
839
841
842 // Set advection velocities
843 StdRegions::VarCoeffType varcoefftypes[] = {
846 for (int j = 0; j < m_spacedim; j++)
847 {
848 varcoeffs[varcoefftypes[j]] = Velocity[j];
849 }
850
851 // Note we are using the
852 // LinearAdvectionDiffusionReaction solver here
853 // instead of HelmSolve since m_filtWidth is negative and
854 // so matrices are not positive definite. Ideally
855 // should allow for negative m_filtWidth coefficient in
856 // HelmSolve
857 // m_recvField->LinearAdvectionDiffusionReactionSolve(
858 // Velocity, forcing, tmpC, -m_filtWidth);
859 m_recvField->LinearAdvectionDiffusionReactionSolve(
860 forcing, tmpC, factors, varcoeffs, varfactors);
861
862 m_evalField->BwdTrans(tmpC, field[i]);
863 timer3.Stop();
864
865 if (m_evalField->GetComm()->GetRank() == 0 &&
866 m_evalField->GetSession()->DefinesCmdLineArgument(
867 "verbose"))
868 {
869 std::cout << "Smoother time (" << m_recvFieldNames[i]
870 << "): " << timer3.TimePerTest(1) << std::endl;
871 }
872 }
873 }
874 else
875 {
876 for (int i = 0; i < m_nRecvVars; ++i)
877 {
878 m_recvField->FwdTrans(rVals[i], tmpC);
879 m_evalField->BwdTrans(tmpC, field[i]);
880 }
881 }
882 timer1.Stop();
883
884 if (m_evalField->GetComm()->GetRank() == 0 &&
885 m_evalField->GetSession()->DefinesCmdLineArgument("verbose"))
886 {
887 std::cout << "Receive total time: " << timer1.TimePerTest(1)
888 << ", ";
889 std::cout << "Receive waiting time: " << timer2.TimePerTest(1)
890 << std::endl;
891 }
892
893 ReceiveStart();
894 }
895}
896
899{
900 LibUtilities::Timer timer1, timer2;
901 timer1.Start();
902
903 int totvars = 3 + m_nRecvVars;
904 int nNotLoc = notLoc.size();
905 int nranks = m_evalField->GetSession()->GetComm()->GetSize();
906
907 Array<OneD, int> thisNLoc(1, m_nPoints - nNotLoc);
908 Array<OneD, int> sizeMap(nranks);
909 m_evalField->GetSession()->GetComm()->AllGather(thisNLoc, sizeMap);
910
911 Array<OneD, int> offsetMap(nranks);
912 offsetMap[0] = 0;
913 for (int i = 1; i < nranks; ++i)
914 {
915 offsetMap[i] = offsetMap[i - 1] + sizeMap[i - 1];
916 }
917 int totNLoc = offsetMap[nranks - 1] + sizeMap[nranks - 1];
918
919 Array<OneD, Array<OneD, NekDouble>> allVals(totvars);
920 for (int i = 0; i < 3; ++i)
921 {
922 allVals[i] = Array<OneD, NekDouble>(m_nPoints);
923 }
924 m_recvField->GetCoords(allVals[0], allVals[1], allVals[2]);
925
926 if (m_spacedim < 3)
927 {
928 Vmath::Zero(m_nPoints, allVals[2], 1);
929 }
930 if (m_spacedim < 2)
931 {
932 Vmath::Zero(m_nPoints, allVals[1], 1);
933 }
934
935 for (int i = 0; i < m_nRecvVars; ++i)
936 {
937 allVals[3 + i] = rVals[i];
938 }
939
940 Array<OneD, Array<OneD, NekDouble>> locatedVals(totvars);
941 for (int i = 0; i < totvars; ++i)
942 {
943 locatedVals[i] = Array<OneD, NekDouble>(totNLoc, -42.0);
944 }
945
946 // only copy points from allVals to locatedVals that were located
947 int offset = offsetMap[m_evalField->GetSession()->GetComm()->GetRank()];
948 for (int i = 0, intPos = 0, locPos = 0; i < m_nPoints; ++i)
949 {
950 if (locPos < nNotLoc && notLoc[locPos] == i)
951 {
952 // do nothing
953 locPos++;
954 }
955 else
956 {
957 for (int k = 0; k < totvars; ++k)
958 {
959 locatedVals[k][offset + intPos] = allVals[k][i];
960 }
961 intPos++;
962 }
963 }
964
965 // send all located points to all ranks. This is probably horribly expensive
966 timer2.Start();
967 for (int i = 0; i < totvars; ++i)
968 {
969 m_evalField->GetSession()->GetComm()->AllGatherv(locatedVals[i],
970 sizeMap, offsetMap);
971 }
972 timer2.Stop();
973
974 if (nNotLoc > 0)
975 {
978 3, locatedVals);
979
980 Array<OneD, Array<OneD, NekDouble>> notLocVals(totvars);
981 for (int j = 0; j < totvars; ++j)
982 {
983 notLocVals[j] = Array<OneD, NekDouble>(nNotLoc);
984 for (int i = 0; i < nNotLoc; ++i)
985 {
986 notLocVals[j][i] = allVals[j][notLoc[i]];
987 }
988 }
991 3, notLocVals);
992
993 // perform a nearest neighbour interpolation from locatedVals to the not
994 // located rVals
996 {
998 std::vector<MultiRegions::ExpListSharedPtr>>>::
999 AllocateSharedPtr(LibUtilities::eNearestNeighbour);
1000 m_extrapInterpolator->CalcWeights(locatedPts, notlocPts);
1001 m_extrapInterpolator->PrintStatistics();
1002 }
1003 m_extrapInterpolator->Interpolate(locatedPts, notlocPts);
1004
1005 for (int j = 3; j < totvars; ++j)
1006 {
1007 for (int i = 0; i < nNotLoc; ++i)
1008 {
1009 allVals[j][notLoc[i]] = notLocVals[j][i];
1010 }
1011 }
1012 }
1013
1014 timer1.Stop();
1015 if (m_evalField->GetComm()->GetRank() == 0 &&
1016 m_evalField->GetSession()->DefinesCmdLineArgument("verbose"))
1017 {
1018 std::cout << "ExtrapolateFields total time: " << timer1.TimePerTest(1);
1019 std::cout << " (AllGathervI: " << timer2.TimePerTest(1) << ")"
1020 << std::endl;
1021 }
1022}
1023
1026{
1027 LibUtilities::Timer timer1;
1028 timer1.Start();
1029
1030#if (defined _WIN32 && _MSC_VER < 1900)
1031 // We need this to make sure boost::format has always
1032 // two digits in the exponents of Scientific notation.
1033 unsigned int old_exponent_format;
1034 old_exponent_format = _set_output_format(_TWO_DIGIT_EXPONENT);
1035 filename = boost::str(boost::format(filename) % m_time);
1036 _set_output_format(old_exponent_format);
1037#else
1038 std::string filename =
1039 boost::str(boost::format(m_config["DUMPRAW"]) % time);
1040#endif
1041
1043 for (int i = 0; i < 3; ++i)
1044 {
1045 tmp[i] = Array<OneD, NekDouble>(m_recvField->GetTotPoints(), 0.0);
1046 }
1047 m_recvField->GetCoords(tmp[0], tmp[1], tmp[2]);
1048
1049 for (int i = 0; i < m_nRecvVars; ++i)
1050 {
1051 tmp[m_spacedim + i] = rVals[i];
1052 }
1053
1054 LibUtilities::CsvIO csvIO(m_evalField->GetSession()->GetComm());
1058 csvIO.Write(filename, rvPts);
1059
1060 timer1.Stop();
1061 if (m_evalField->GetComm()->GetRank() == 0 &&
1062 m_evalField->GetSession()->DefinesCmdLineArgument("verbose"))
1063 {
1064 std::cout << "DumpRawFields total time: " << timer1.TimePerTest(1)
1065 << std::endl;
1066 }
1067}
1068} // namespace Nektar::SolverUtils
#define OUTPUT_FREQ
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:242
A class that contains algorithms for interpolation between pts fields, expansions and different meshe...
void Write(const std::string &outFile, const PtsFieldSharedPtr &ptsField, const bool backup=false)
Save a pts field to a file.
Definition: CsvIO.cpp:68
NekDouble TimePerTest(unsigned int n)
Returns amount of seconds per iteration in a test with n iterations.
Definition: Timer.cpp:65
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
static SOLVER_UTILS_EXPORT void InterpCallback(const int entities_dim, const int n_local_vertex, const int n_local_element, const int n_local_polhyedra, const int n_distant_point, const double local_coordinates[], const int local_connectivity_index[], const int local_connectivity[], const int local_polyhedra_face_index[], const int local_polyhedra_cell_to_face_connectivity[], const int local_polyhedra_face_connectivity_index[], const int local_polyhedra_face_connectivity[], const double distant_points_coordinates[], const int distant_points_location[], const float distant_points_distance[], const int distant_points_barycentric_coordinates_index[], const double distant_points_barycentric_coordinates[], const int stride, const cwipi_solver_type_t solver_type, const void *local_field, void *distant_field)
Array< OneD, Array< OneD, NekDouble > > m_oldFields
Definition: CouplingCwipi.h:94
SOLVER_UTILS_EXPORT void v_Receive(const int step, const NekDouble time, Array< OneD, Array< OneD, NekDouble > > &field, std::vector< std::string > &varNames) override
void AddElementsToMesh(T geom, int &coordsPos, int &connecPos, int &conidxPos)
Array< OneD, Array< OneD, NekDouble > > m_sendField
Definition: CouplingCwipi.h:90
SOLVER_UTILS_EXPORT CouplingCwipi(MultiRegions::ExpListSharedPtr field)
SOLVER_UTILS_EXPORT void v_Init() override
void ExtrapolateFields(Array< OneD, Array< OneD, NekDouble > > &rVals, Array< OneD, int > &notLoc)
std::shared_ptr< FieldUtils::Interpolator< std::vector< MultiRegions::ExpListSharedPtr > > > m_sendInterpolator
void EvaluateFields(Array< OneD, Array< OneD, NekDouble > > interpField, Array< OneD, Array< OneD, NekDouble > > distCoords)
SOLVER_UTILS_EXPORT void v_Finalize() override
SOLVER_UTILS_EXPORT ~CouplingCwipi() override
static CouplingSharedPtr create(MultiRegions::ExpListSharedPtr field)
Creates an instance of this class.
Definition: CouplingCwipi.h:58
Array< OneD, Array< OneD, NekDouble > > m_newFields
Definition: CouplingCwipi.h:95
SOLVER_UTILS_EXPORT void v_Send(const int step, const NekDouble time, const Array< OneD, const Array< OneD, NekDouble > > &field, std::vector< std::string > &varNames) override
SOLVER_UTILS_EXPORT void SendCallback(Array< OneD, Array< OneD, NekDouble > > &interpField, Array< OneD, Array< OneD, NekDouble > > &distCoords)
void ReadConfig(LibUtilities::SessionReaderSharedPtr session)
void ReceiveCwipi(const int step, const NekDouble time, Array< OneD, Array< OneD, NekDouble > > &field)
std::shared_ptr< FieldUtils::Interpolator< std::vector< MultiRegions::ExpListSharedPtr > > > m_extrapInterpolator
void DumpRawFields(const NekDouble time, Array< OneD, Array< OneD, NekDouble > > rVals)
MultiRegions::ExpListSharedPtr m_recvField
Definition: CouplingCwipi.h:92
virtual SOLVER_UTILS_EXPORT void v_Init()
Definition: Coupling.cpp:57
MultiRegions::ExpListSharedPtr m_evalField
Definition: Coupling.h:107
std::vector< std::string > m_sendFieldNames
Definition: Coupling.h:110
CouplingConfigMap m_config
Definition: Coupling.h:105
std::vector< std::string > m_recvFieldNames
Definition: Coupling.h:114
SOLVER_UTILS_EXPORT std::vector< int > GenerateVariableMapping(std::vector< std::string > &vars, std::vector< std::string > &transVars)
Definition: Coupling.cpp:131
static MeshGraphSharedPtr Read(const LibUtilities::SessionReaderSharedPtr pSession, LibUtilities::DomainRangeShPtr rng=LibUtilities::NullDomainRangeShPtr, bool fillGraph=true, SpatialDomains::MeshGraphSharedPtr partitionedGraph=nullptr)
Definition: MeshGraphIO.cpp:51
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< PtsField > PtsFieldSharedPtr
Definition: PtsField.h:184
static VarFactorsMap NullVarFactorsMap
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
std::map< StdRegions::ConstFactorType, Array< OneD, NekDouble > > VarFactorsMap
static const NekDouble kNekZeroTol
static std::map< std::string, SendCallbackType > SendCallbackMap
CouplingFactory & GetCouplingFactory()
Declaration of the Coupling factory singleton.
Definition: Coupling.cpp:40
std::map< int, TriGeomSharedPtr > TriGeomMap
Definition: TriGeom.h:57
std::map< int, PyrGeomSharedPtr > PyrGeomMap
Definition: PyrGeom.h:76
std::map< int, QuadGeomSharedPtr > QuadGeomMap
Definition: QuadGeom.h:53
std::map< int, SegGeomSharedPtr > SegGeomMap
Definition: SegGeom.h:50
std::map< int, TetGeomSharedPtr > TetGeomMap
Definition: TetGeom.h:84
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:174
std::map< int, PrismGeomSharedPtr > PrismGeomMap
Definition: PrismGeom.h:83
std::shared_ptr< PointGeom > PointGeomSharedPtr
Definition: Geometry.h:57
std::map< int, HexGeomSharedPtr > HexGeomMap
Definition: HexGeom.h:85
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:430
std::map< StdRegions::VarCoeffType, VarCoeffEntry > VarCoeffMap
Definition: StdRegions.hpp:375
StdRegions::ConstFactorMap factors
double NekDouble
void Svtsvtp(int n, const T alpha, const T *x, int incx, const T beta, const T *y, int incy, T *z, int incz)
Svtsvtp (scalar times vector plus scalar times vector):
Definition: Vmath.hpp:473
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
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.hpp:273
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825