46#include <boost/format.hpp>
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)
79 for (
int i = 0; i < n_distant_point; ++i)
82 for (
int j = 0; j < 3; ++j)
84 distCoords[i][j] = distant_points_coordinates[3 * i + j];
88 std::stringstream sst;
89 sst << entities_dim <<
"," << n_local_vertex <<
"," << stride;
92 ASSERTL0(interpField.size() == stride,
"size mismatch");
93 ASSERTL0(interpField[0].size() == n_distant_point,
"size mismatch");
95 for (
int i = 0; i < n_distant_point; i++)
97 for (
int j = 0; j < stride; ++j)
99 ((
double *)distant_field)[i * stride + j] = interpField[j][i];
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)
117 m_config[
"SENDMETHOD"] =
"EVALUATE";
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",
142 cwipi_add_local_int_control_parameter(
"receiveTag",
m_recvTag);
153 cwipi_solver_type_t solver_type = CWIPI_SOLVER_CELL_VERTEX;
155 cwipi_create_coupling(
156 m_couplingName.c_str(), CWIPI_COUPLING_PARALLEL_WITH_PARTITIONING,
159 cwipi_synchronize_control_parameter(
m_config[
"REMOTENAME"].c_str());
162 m_evalField->GetSession()->DefinesCmdLineArgument(
"verbose"))
164 cwipi_dump_application_properties();
167 m_sendTag = cwipi_get_distant_int_control_parameter(
168 m_config[
"REMOTENAME"].c_str(),
"receiveTag");
170 if (cwipi_has_int_parameter(
m_config[
"REMOTENAME"].c_str(),
"nRecvVars"))
172 int remoteNRecvVars = cwipi_get_distant_int_control_parameter(
173 m_config[
"REMOTENAME"].c_str(),
"nRecvVars");
175 "Number of local send vars different to remote received vars");
178 if (cwipi_has_int_parameter(
m_config[
"REMOTENAME"].c_str(),
"nSendVars"))
180 int remoteNSendVars = cwipi_get_distant_int_control_parameter(
181 m_config[
"REMOTENAME"].c_str(),
"nSendVars");
183 "Number of local receive vars different to remote sent vars");
194 m_evalField->GetSession()->DefinesCmdLineArgument(
"verbose"))
196 std::cout <<
"locating..." << std::endl;
223 ASSERTL0(session->DefinesElement(
"Nektar/Coupling"),
224 "No Coupling config found");
226 m_config[
"LOCALNAME"] = session->GetCmdLineArgument<std::string>(
"cwipi");
228 TiXmlElement *vCoupling = session->GetElement(
"Nektar/Coupling");
229 ASSERTL0(vCoupling,
"Invalid Coupling config");
236 int oversamp = std::stoi(
m_config[
"OVERSAMPLE"]);
240 recvGraph->SetExpansionInfoToPointOrder(
241 oversamp +
m_evalField->GetExp(0)->GetNumPoints(0));
245 m_evalField->GetSession(), recvGraph,
"DefaultVar");
261 m_recvField->GetCoords(coords[0], coords[1], coords[2]);
268 m_points[3 * i + 0] = double(coords[0][i]);
272 m_points[3 * i + 1] = double(coords[1][i]);
281 m_points[3 * i + 2] = double(coords[2][i]);
308 std::stringstream sst;
313 std::placeholders::_2);
322 int nOutPts = distCoords.size();
325 for (
int i = 0; i < nOutPts; ++i)
332 elmtid =
m_evalField->GetExpIndex(distCoords[i], Lcoords, tol);
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]) +
")");
363 ASSERTL0(!(boost::math::isnan)(value),
"new value is not a number");
364 interpField[f][i] = value;
371 const double *distCoords =
376 for (
int i = 0; i < 3; ++i)
380 m_evalField->GetCoords(local[0], local[1], local[2]);
385 for (
int i = 0; i < 3; ++i)
388 for (
int j = 0; j < nPts; ++j)
390 dist[i][j] = distCoords[3 * j + i];
397 if (boost::to_upper_copy(
m_config[
"SENDMETHOD"]) ==
"SHEPARD")
421 seggeom = graph->GetAllSegGeoms();
425 trigeom = graph->GetAllTriGeoms();
426 quadgeom = graph->GetAllQuadGeoms();
430 tetgeom = graph->GetAllTetGeoms();
431 pyrgeom = graph->GetAllPyrGeoms();
432 prismgeom = graph->GetAllPrismGeoms();
433 hexgeom = graph->GetAllHexGeoms();
436 int nVerts = graph->GetNvertices();
437 int nElts = seggeom.size() + trigeom.size() + quadgeom.size() +
438 tetgeom.size() + pyrgeom.size() + prismgeom.size() +
442 m_coords = (
double *)malloc(
sizeof(
double) * 3 * nVerts);
444 int tmp = 2 * seggeom.size() + 3 * trigeom.size() + 4 * quadgeom.size() +
445 4 * tetgeom.size() + 5 * pyrgeom.size() + 6 * prismgeom.size() +
447 m_connec = (
int *)malloc(
sizeof(
int) * tmp);
449 m_connecIdx = (
int *)malloc(
sizeof(
int) * (nElts + 1));
502 int kNverts = T::mapped_type::element_type::kNverts;
505 for (
auto it = geom.begin(); it != geom.end(); it++)
508 for (
int j = 0; j < kNverts; ++j)
510 vert = it->second->GetVertex(j);
511 vertID = vert->GetVid();
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]);
547 if (boost::to_upper_copy(
m_config[
"SENDMETHOD"]) ==
"NEARESTNEIGHBOUR" ||
548 boost::to_upper_copy(
m_config[
"SENDMETHOD"]) ==
"SHEPARD")
572 std::vector<std::string> &varNames)
585 std::vector<int> sendVarsToVars =
588 for (
int i = 0; i < sendVarsToVars.size(); ++i)
594 m_evalField->GetSession()->DefinesCmdLineArgument(
"verbose"))
596 std::cout <<
"sending fields at i = " << step <<
", t = " << time
604 strcpy(sendFN,
"dummyName");
611 m_evalField->GetSession()->DefinesCmdLineArgument(
"verbose"))
613 std::cout <<
"Send total time: " << timer1.
TimePerTest(1)
635 m_evalField->GetSession()->DefinesCmdLineArgument(
"verbose"))
637 std::cout <<
"Send waiting time: " << timer1.
TimePerTest(1)
654 strcpy(recFN,
"dummyName");
661 m_evalField->GetSession()->DefinesCmdLineArgument(
"verbose"))
663 std::cout <<
"Receive start time: " << timer1.
TimePerTest(1)
670 std::vector<std::string> &varNames)
678 std::vector<int> recvVarsToVars =
681 for (
int i = 0; i < recvVarsToVars.size(); ++i)
683 recvFields[i] =
field[recvVarsToVars[i]];
714 1, recvFields[i], 1);
733 m_evalField->GetSession()->DefinesCmdLineArgument(
"verbose"))
735 std::cout <<
"waiting for receive at i = " << step
736 <<
", t = " << time << std::endl;
756 int nNotLoc = cwipi_get_n_not_located_points(
m_couplingName.c_str());
760 std::cout <<
"WARNING: relocating " << nNotLoc <<
" of "
766 for (
int i = 0; i < nNotLoc; ++i)
768 notLoc[i] = tmp[i] - 1;
771 if (boost::to_upper_copy(
m_config[
"NOTLOCMETHOD"]) ==
"KEEP")
782 for (
int i = 0, locPos = 0, intPos = 0; i <
m_nPoints; ++i)
784 if (locPos < nNotLoc && notLoc[locPos] == i)
799 if (boost::to_upper_copy(
m_config[
"NOTLOCMETHOD"]) ==
"EXTRAPOLATE")
801 int doExtrapolate = 0;
808 if (doExtrapolate > 0)
848 varcoeffs[varcoefftypes[j]] = Velocity[j];
859 m_recvField->LinearAdvectionDiffusionReactionSolve(
860 forcing, tmpC,
factors, varcoeffs, varfactors);
885 m_evalField->GetSession()->DefinesCmdLineArgument(
"verbose"))
887 std::cout <<
"Receive total time: " << timer1.
TimePerTest(1)
889 std::cout <<
"Receive waiting time: " << timer2.
TimePerTest(1)
904 int nNotLoc = notLoc.size();
905 int nranks =
m_evalField->GetSession()->GetComm()->GetSize();
909 m_evalField->GetSession()->GetComm()->AllGather(thisNLoc, sizeMap);
913 for (
int i = 1; i < nranks; ++i)
915 offsetMap[i] = offsetMap[i - 1] + sizeMap[i - 1];
917 int totNLoc = offsetMap[nranks - 1] + sizeMap[nranks - 1];
920 for (
int i = 0; i < 3; ++i)
924 m_recvField->GetCoords(allVals[0], allVals[1], allVals[2]);
937 allVals[3 + i] = rVals[i];
941 for (
int i = 0; i < totvars; ++i)
947 int offset = offsetMap[
m_evalField->GetSession()->GetComm()->GetRank()];
948 for (
int i = 0, intPos = 0, locPos = 0; i <
m_nPoints; ++i)
950 if (locPos < nNotLoc && notLoc[locPos] == i)
957 for (
int k = 0; k < totvars; ++k)
959 locatedVals[k][offset + intPos] = allVals[k][i];
967 for (
int i = 0; i < totvars; ++i)
969 m_evalField->GetSession()->GetComm()->AllGatherv(locatedVals[i],
981 for (
int j = 0; j < totvars; ++j)
984 for (
int i = 0; i < nNotLoc; ++i)
986 notLocVals[j][i] = allVals[j][notLoc[i]];
998 std::vector<MultiRegions::ExpListSharedPtr>>>::
1005 for (
int j = 3; j < totvars; ++j)
1007 for (
int i = 0; i < nNotLoc; ++i)
1009 allVals[j][notLoc[i]] = notLocVals[j][i];
1016 m_evalField->GetSession()->DefinesCmdLineArgument(
"verbose"))
1018 std::cout <<
"ExtrapolateFields total time: " << timer1.
TimePerTest(1);
1019 std::cout <<
" (AllGathervI: " << timer2.
TimePerTest(1) <<
")"
1030#if (defined _WIN32 && _MSC_VER < 1900)
1033 unsigned int old_exponent_format;
1034 old_exponent_format = _set_output_format(_TWO_DIGIT_EXPONENT);
1036 _set_output_format(old_exponent_format);
1038 std::string filename =
1043 for (
int i = 0; i < 3; ++i)
1058 csvIO.
Write(filename, rvPts);
1062 m_evalField->GetSession()->DefinesCmdLineArgument(
"verbose"))
1064 std::cout <<
"DumpRawFields total time: " << timer1.
TimePerTest(1)
#define ASSERTL0(condition, msg)
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
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.
NekDouble TimePerTest(unsigned int n)
Returns amount of seconds per iteration in a test with n iterations.
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)
static std::string className
Array< OneD, Array< OneD, NekDouble > > m_oldFields
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
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 > ¬Loc)
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.
Array< OneD, Array< OneD, NekDouble > > m_newFields
void SetupSendInterpolation()
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
std::map< int, int > m_vertMap
void DumpRawFields(const NekDouble time, Array< OneD, Array< OneD, NekDouble > > rVals)
MultiRegions::ExpListSharedPtr m_recvField
virtual SOLVER_UTILS_EXPORT void v_Init()
MultiRegions::ExpListSharedPtr m_evalField
std::string m_couplingName
std::vector< std::string > m_sendFieldNames
CouplingConfigMap m_config
std::vector< std::string > m_recvFieldNames
SOLVER_UTILS_EXPORT std::vector< int > GenerateVariableMapping(std::vector< std::string > &vars, std::vector< std::string > &transVars)
static MeshGraphSharedPtr Read(const LibUtilities::SessionReaderSharedPtr pSession, LibUtilities::DomainRangeShPtr rng=LibUtilities::NullDomainRangeShPtr, bool fillGraph=true, SpatialDomains::MeshGraphSharedPtr partitionedGraph=nullptr)
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< PtsField > PtsFieldSharedPtr
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.
std::map< int, TriGeomSharedPtr > TriGeomMap
std::map< int, PyrGeomSharedPtr > PyrGeomMap
std::map< int, QuadGeomSharedPtr > QuadGeomMap
std::map< int, SegGeomSharedPtr > SegGeomMap
std::map< int, TetGeomSharedPtr > TetGeomMap
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
std::map< int, PrismGeomSharedPtr > PrismGeomMap
std::shared_ptr< PointGeom > PointGeomSharedPtr
std::map< int, HexGeomSharedPtr > HexGeomMap
std::map< ConstFactorType, NekDouble > ConstFactorMap
std::map< StdRegions::VarCoeffType, VarCoeffEntry > VarCoeffMap
StdRegions::ConstFactorMap factors
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):
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
void Zero(int n, T *x, const int incx)
Zero vector.
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)