45#include <boost/format.hpp>
46#include <boost/functional/hash.hpp>
60 const int entities_dim,
const int n_local_vertex,
61 [[maybe_unused]]
const int n_local_element,
62 [[maybe_unused]]
const int n_local_polyhedra,
const int n_distant_point,
63 [[maybe_unused]]
const double local_coordinates[],
64 [[maybe_unused]]
const int local_connectivity_index[],
65 [[maybe_unused]]
const int local_connectivity[],
66 [[maybe_unused]]
const int local_polyhedra_face_index[],
67 [[maybe_unused]]
const int local_polyhedra_cell_to_face_connectivity[],
68 [[maybe_unused]]
const int local_polyhedra_face_connectivity_index[],
69 [[maybe_unused]]
const int local_polyhedra_face_connectivity[],
70 const double distant_points_coordinates[],
71 [[maybe_unused]]
const int distant_points_location[],
72 [[maybe_unused]]
const float distant_points_distance[],
73 [[maybe_unused]]
const int distant_points_barycentric_coordinates_index[],
74 [[maybe_unused]]
const double distant_points_barycentric_coordinates[],
75 const int stride, [[maybe_unused]]
const cwipi_solver_type_t solver_type,
76 [[maybe_unused]]
const void *local_field,
void *distant_field)
81 for (
int i = 0; i < n_distant_point; ++i)
84 for (
int j = 0; j < 3; ++j)
86 distCoords[i][j] = distant_points_coordinates[3 * i + j];
90 std::stringstream sst;
91 sst << entities_dim <<
"," << n_local_vertex <<
"," << stride;
94 ASSERTL0(interpField.size() == stride,
"size mismatch");
95 ASSERTL0(interpField[0].size() == n_distant_point,
"size mismatch");
97 for (
int i = 0; i < n_distant_point; i++)
99 for (
int j = 0; j < stride; ++j)
101 ((
double *)distant_field)[i * stride + j] = interpField[j][i];
107 :
Coupling(
field), m_sendHandle(-1), m_recvHandle(-1), m_lastSend(-1E6),
108 m_lastReceive(-1E6), m_points(nullptr), m_coords(nullptr),
109 m_connecIdx(nullptr), m_connec(nullptr), m_rValsInterl(nullptr),
110 m_sValsInterl(nullptr)
119 m_config[
"SENDMETHOD"] =
"EVALUATE";
129 cwipi_add_local_int_control_parameter(
"nSendVars",
m_nSendVars);
130 cwipi_add_local_int_control_parameter(
"nRecvVars",
m_nRecvVars);
131 cwipi_add_local_string_control_parameter(
132 "recvFieldNames",
m_config[
"RECEIVEVARIABLES"].c_str());
133 cwipi_add_local_string_control_parameter(
"sendFieldNames",
144 cwipi_add_local_int_control_parameter(
"receiveTag",
m_recvTag);
155 cwipi_solver_type_t solver_type = CWIPI_SOLVER_CELL_VERTEX;
157 cwipi_create_coupling(
158 m_couplingName.c_str(), CWIPI_COUPLING_PARALLEL_WITH_PARTITIONING,
161 cwipi_synchronize_control_parameter(
m_config[
"REMOTENAME"].c_str());
164 m_evalField->GetSession()->DefinesCmdLineArgument(
"verbose"))
166 cwipi_dump_application_properties();
169 m_sendTag = cwipi_get_distant_int_control_parameter(
170 m_config[
"REMOTENAME"].c_str(),
"receiveTag");
172 if (cwipi_has_int_parameter(
m_config[
"REMOTENAME"].c_str(),
"nRecvVars"))
174 int remoteNRecvVars = cwipi_get_distant_int_control_parameter(
175 m_config[
"REMOTENAME"].c_str(),
"nRecvVars");
177 "Number of local send vars different to remote received vars");
180 if (cwipi_has_int_parameter(
m_config[
"REMOTENAME"].c_str(),
"nSendVars"))
182 int remoteNSendVars = cwipi_get_distant_int_control_parameter(
183 m_config[
"REMOTENAME"].c_str(),
"nSendVars");
185 "Number of local receive vars different to remote sent vars");
196 m_evalField->GetSession()->DefinesCmdLineArgument(
"verbose"))
198 cout <<
"locating..." << endl;
225 ASSERTL0(session->DefinesElement(
"Nektar/Coupling"),
226 "No Coupling config found");
228 m_config[
"LOCALNAME"] = session->GetCmdLineArgument<std::string>(
"cwipi");
230 TiXmlElement *vCoupling = session->GetElement(
"Nektar/Coupling");
231 ASSERTL0(vCoupling,
"Invalid Coupling config");
238 int oversamp = boost::lexical_cast<int>(
m_config[
"OVERSAMPLE"]);
242 recvGraph->SetExpansionInfoToPointOrder(
243 oversamp +
m_evalField->GetExp(0)->GetNumPoints(0));
247 m_evalField->GetSession(), recvGraph,
"DefaultVar");
263 m_recvField->GetCoords(coords[0], coords[1], coords[2]);
270 m_points[3 * i + 0] = double(coords[0][i]);
274 m_points[3 * i + 1] = double(coords[1][i]);
283 m_points[3 * i + 2] = double(coords[2][i]);
310 std::stringstream sst;
315 std::placeholders::_2);
324 int nOutPts = distCoords.size();
327 for (
int i = 0; i < nOutPts; ++i)
334 elmtid =
m_evalField->GetExpIndex(distCoords[i], Lcoords, tol);
353 "no element found for (" +
354 boost::lexical_cast<string>(distCoords[i][0]) +
", " +
355 boost::lexical_cast<string>(distCoords[i][1]) +
", " +
356 boost::lexical_cast<string>(distCoords[i][2]) +
")");
365 ASSERTL0(!(boost::math::isnan)(value),
"new value is not a number");
366 interpField[f][i] = value;
373 const double *distCoords =
378 for (
int i = 0; i < 3; ++i)
382 m_evalField->GetCoords(local[0], local[1], local[2]);
387 for (
int i = 0; i < 3; ++i)
390 for (
int j = 0; j < nPts; ++j)
392 dist[i][j] = distCoords[3 * j + i];
399 if (boost::to_upper_copy(
m_config[
"SENDMETHOD"]) ==
"SHEPARD")
423 seggeom = graph->GetAllSegGeoms();
427 trigeom = graph->GetAllTriGeoms();
428 quadgeom = graph->GetAllQuadGeoms();
432 tetgeom = graph->GetAllTetGeoms();
433 pyrgeom = graph->GetAllPyrGeoms();
434 prismgeom = graph->GetAllPrismGeoms();
435 hexgeom = graph->GetAllHexGeoms();
438 int nVerts = graph->GetNvertices();
439 int nElts = seggeom.size() + trigeom.size() + quadgeom.size() +
440 tetgeom.size() + pyrgeom.size() + prismgeom.size() +
444 m_coords = (
double *)malloc(
sizeof(
double) * 3 * nVerts);
446 int tmp = 2 * seggeom.size() + 3 * trigeom.size() + 4 * quadgeom.size() +
447 4 * tetgeom.size() + 5 * pyrgeom.size() + 6 * prismgeom.size() +
449 m_connec = (
int *)malloc(
sizeof(
int) * tmp);
451 m_connecIdx = (
int *)malloc(
sizeof(
int) * (nElts + 1));
503 int kNverts = T::mapped_type::element_type::kNverts;
506 for (
auto it = geom.begin(); it != geom.end(); it++)
509 for (
int j = 0; j < kNverts; ++j)
511 vert = it->second->GetVertex(j);
512 vertID = vert->GetVid();
518 vert->GetCoords(x[0], x[1], x[2]);
519 m_coords[3 * coordsPos + 0] = double(x[0]);
520 m_coords[3 * coordsPos + 1] = double(x[1]);
521 m_coords[3 * coordsPos + 2] = double(x[2]);
548 if (boost::to_upper_copy(
m_config[
"SENDMETHOD"]) ==
"NEARESTNEIGHBOUR" ||
549 boost::to_upper_copy(
m_config[
"SENDMETHOD"]) ==
"SHEPARD")
573 vector<string> &varNames)
586 vector<int> sendVarsToVars =
589 for (
int i = 0; i < sendVarsToVars.size(); ++i)
595 m_evalField->GetSession()->DefinesCmdLineArgument(
"verbose"))
597 cout <<
"sending fields at i = " << step <<
", t = " << time
605 strcpy(sendFN,
"dummyName");
612 m_evalField->GetSession()->DefinesCmdLineArgument(
"verbose"))
614 cout <<
"Send total time: " << timer1.
TimePerTest(1) << endl;
635 m_evalField->GetSession()->DefinesCmdLineArgument(
"verbose"))
637 cout <<
"Send waiting time: " << timer1.
TimePerTest(1) << endl;
653 strcpy(recFN,
"dummyName");
660 m_evalField->GetSession()->DefinesCmdLineArgument(
"verbose"))
662 cout <<
"Receive start time: " << timer1.
TimePerTest(1) << endl;
668 vector<string> &varNames)
676 vector<int> recvVarsToVars =
679 for (
int i = 0; i < recvVarsToVars.size(); ++i)
681 recvFields[i] =
field[recvVarsToVars[i]];
712 1, recvFields[i], 1);
731 m_evalField->GetSession()->DefinesCmdLineArgument(
"verbose"))
733 cout <<
"waiting for receive at i = " << step <<
", t = " << time
754 int nNotLoc = cwipi_get_n_not_located_points(
m_couplingName.c_str());
758 cout <<
"WARNING: relocating " << nNotLoc <<
" of " <<
m_nPoints
759 <<
" points" << endl;
764 for (
int i = 0; i < nNotLoc; ++i)
766 notLoc[i] = tmp[i] - 1;
769 if (boost::to_upper_copy(
m_config[
"NOTLOCMETHOD"]) ==
"KEEP")
780 for (
int i = 0, locPos = 0, intPos = 0; i <
m_nPoints; ++i)
782 if (locPos < nNotLoc && notLoc[locPos] == i)
797 if (boost::to_upper_copy(
m_config[
"NOTLOCMETHOD"]) ==
"EXTRAPOLATE")
799 int doExtrapolate = 0;
806 if (doExtrapolate > 0)
846 varcoeffs[varcoefftypes[j]] = Velocity[j];
857 m_recvField->LinearAdvectionDiffusionReactionSolve(
858 forcing, tmpC,
factors, varcoeffs, varfactors);
883 m_evalField->GetSession()->DefinesCmdLineArgument(
"verbose"))
885 cout <<
"Receive total time: " << timer1.
TimePerTest(1) <<
", ";
886 cout <<
"Receive waiting time: " << timer2.
TimePerTest(1) << endl;
900 int nNotLoc = notLoc.size();
901 int nranks =
m_evalField->GetSession()->GetComm()->GetSize();
905 m_evalField->GetSession()->GetComm()->AllGather(thisNLoc, sizeMap);
909 for (
int i = 1; i < nranks; ++i)
911 offsetMap[i] = offsetMap[i - 1] + sizeMap[i - 1];
913 int totNLoc = offsetMap[nranks - 1] + sizeMap[nranks - 1];
916 for (
int i = 0; i < 3; ++i)
920 m_recvField->GetCoords(allVals[0], allVals[1], allVals[2]);
933 allVals[3 + i] = rVals[i];
937 for (
int i = 0; i < totvars; ++i)
943 int offset = offsetMap[
m_evalField->GetSession()->GetComm()->GetRank()];
944 for (
int i = 0, intPos = 0, locPos = 0; i <
m_nPoints; ++i)
946 if (locPos < nNotLoc && notLoc[locPos] == i)
953 for (
int k = 0; k < totvars; ++k)
955 locatedVals[k][offset + intPos] = allVals[k][i];
963 for (
int i = 0; i < totvars; ++i)
965 m_evalField->GetSession()->GetComm()->AllGatherv(locatedVals[i],
977 for (
int j = 0; j < totvars; ++j)
980 for (
int i = 0; i < nNotLoc; ++i)
982 notLocVals[j][i] = allVals[j][notLoc[i]];
994 std::vector<MultiRegions::ExpListSharedPtr>>>::
1001 for (
int j = 3; j < totvars; ++j)
1003 for (
int i = 0; i < nNotLoc; ++i)
1005 allVals[j][notLoc[i]] = notLocVals[j][i];
1012 m_evalField->GetSession()->DefinesCmdLineArgument(
"verbose"))
1014 cout <<
"ExtrapolateFields total time: " << timer1.
TimePerTest(1);
1015 cout <<
" (AllGathervI: " << timer2.
TimePerTest(1) <<
")" << endl;
1025#if (defined _WIN32 && _MSC_VER < 1900)
1028 unsigned int old_exponent_format;
1029 old_exponent_format = _set_output_format(_TWO_DIGIT_EXPONENT);
1031 _set_output_format(old_exponent_format);
1033 std::string filename =
1038 for (
int i = 0; i < 3; ++i)
1053 csvIO.
Write(filename, rvPts);
1057 m_evalField->GetSession()->DefinesCmdLineArgument(
"verbose"))
1059 cout <<
"DumpRawFields total time: " << timer1.
TimePerTest(1) << endl;
#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)