preCICE v3.1.1
Loading...
Searching...
No Matches
RadialBasisFctMapping.hpp
Go to the documentation of this file.
1#pragma once
2
3#include <Eigen/Cholesky>
4#include <Eigen/Core>
5
7#include "com/Extra.hpp"
10#include "mesh/Filter.hpp"
12#include "profiling/Event.hpp"
13#include "utils/IntraComm.hpp"
14
15namespace precice {
16namespace mapping {
17
29template <typename SOLVER_T, typename... Args>
30class RadialBasisFctMapping : public RadialBasisFctBaseMapping<typename SOLVER_T::BASIS_FUNCTION_T> {
31public:
32 using RADIAL_BASIS_FUNCTION_T = typename SOLVER_T::BASIS_FUNCTION_T;
42 Mapping::Constraint constraint,
43 int dimensions,
45 std::array<bool, 3> deadAxis,
46 Polynomial polynomial,
47 Args... args);
48
50 void computeMapping() final override;
51
53 void clear() final override;
54
56 std::string getName() const final override;
57
58private:
59 precice::logging::Logger _log{"mapping::RadialBasisFctMapping"};
60
61 // The actual solver
63
65 void mapConservative(const time::Sample &inData, Eigen::VectorXd &outData) final override;
66
68 void mapConsistent(const time::Sample &inData, Eigen::VectorXd &outData) final override;
69
72
75};
76
77// --------------------------------------------------- HEADER IMPLEMENTATIONS
78
79template <typename SOLVER_T, typename... Args>
81 Mapping::Constraint constraint,
82 int dimensions,
84 std::array<bool, 3> deadAxis,
85 Polynomial polynomial,
86 Args... args)
87 : RadialBasisFctBaseMapping<RADIAL_BASIS_FUNCTION_T>(constraint, dimensions, function, deadAxis, Mapping::InitialGuessRequirement::None),
88 _polynomial(polynomial),
89 optionalArgs(std::make_tuple(std::forward<Args>(args)...))
90{
91 PRECICE_CHECK(!(RADIAL_BASIS_FUNCTION_T::isStrictlyPositiveDefinite() && polynomial == Polynomial::ON), "The integrated polynomial (polynomial=\"on\") is not supported for the selected radial-basis function. Please select another radial-basis function or change the polynomial configuration.");
92}
93
94template <typename SOLVER_T, typename... Args>
96{
98
99 precice::profiling::Event e("map.rbf.computeMapping.From" + this->input()->getName() + "To" + this->output()->getName(), profiling::Synchronize);
100
101 PRECICE_ASSERT(this->input()->getDimensions() == this->output()->getDimensions(),
102 this->input()->getDimensions(), this->output()->getDimensions());
103 PRECICE_ASSERT(this->getDimensions() == this->output()->getDimensions(),
104 this->getDimensions(), this->output()->getDimensions());
105
106 mesh::PtrMesh inMesh;
107 mesh::PtrMesh outMesh;
108
109 if (this->hasConstraint(Mapping::CONSERVATIVE)) {
110 inMesh = this->output();
111 outMesh = this->input();
112 } else { // Consistent or scaled consistent
113 inMesh = this->input();
114 outMesh = this->output();
115 }
116
118
119 // Input mesh may have overlaps
120 mesh::Mesh filteredInMesh("filteredInMesh", inMesh->getDimensions(), mesh::Mesh::MESH_ID_UNDEFINED);
121 mesh::filterMesh(filteredInMesh, *inMesh, [&](const mesh::Vertex &v) { return v.isOwner(); });
122
123 // Send the mesh
126
127 } else { // Parallel Primary rank or Serial
128
129 mesh::Mesh globalInMesh(inMesh->getName(), inMesh->getDimensions(), mesh::Mesh::MESH_ID_UNDEFINED);
130 mesh::Mesh globalOutMesh(outMesh->getName(), outMesh->getDimensions(), mesh::Mesh::MESH_ID_UNDEFINED);
131
133 {
134 // Input mesh may have overlaps
135 mesh::Mesh filteredInMesh("filteredInMesh", inMesh->getDimensions(), mesh::Mesh::MESH_ID_UNDEFINED);
136 mesh::filterMesh(filteredInMesh, *inMesh, [&](const mesh::Vertex &v) { return v.isOwner(); });
137 globalInMesh.addMesh(filteredInMesh);
138 globalOutMesh.addMesh(*outMesh);
139 }
140
141 // Receive mesh
142 for (Rank secondaryRank : utils::IntraComm::allSecondaryRanks()) {
143 mesh::Mesh secondaryInMesh(inMesh->getName(), inMesh->getDimensions(), mesh::Mesh::MESH_ID_UNDEFINED);
144 com::receiveMesh(*utils::IntraComm::getCommunication(), secondaryRank, secondaryInMesh);
145 globalInMesh.addMesh(secondaryInMesh);
146
147 mesh::Mesh secondaryOutMesh(outMesh->getName(), outMesh->getDimensions(), mesh::Mesh::MESH_ID_UNDEFINED);
148 com::receiveMesh(*utils::IntraComm::getCommunication(), secondaryRank, secondaryOutMesh);
149 globalOutMesh.addMesh(secondaryOutMesh);
150 }
151
152 } else { // Serial
153 globalInMesh.addMesh(*inMesh);
154 globalOutMesh.addMesh(*outMesh);
155 }
156
157 // Forwarding the tuples here requires some template magic I don't want to implement
158 if constexpr (std::tuple_size_v<std::tuple<Args...>>> 0) {
159 _rbfSolver = std::make_unique<SOLVER_T>(this->_basisFunction, globalInMesh, boost::irange<Eigen::Index>(0, globalInMesh.nVertices()),
160 globalOutMesh, boost::irange<Eigen::Index>(0, globalOutMesh.nVertices()), this->_deadAxis, _polynomial, std::get<0>(optionalArgs));
161 } else {
162 _rbfSolver = std::make_unique<SOLVER_T>(this->_basisFunction, globalInMesh, boost::irange<Eigen::Index>(0, globalInMesh.nVertices()),
163 globalOutMesh, boost::irange<Eigen::Index>(0, globalOutMesh.nVertices()), this->_deadAxis, _polynomial);
164 }
165 }
166 this->_hasComputedMapping = true;
167 PRECICE_DEBUG("Compute Mapping is Completed.");
168}
169
170template <typename SOLVER_T, typename... Args>
172{
174 _rbfSolver.reset();
175 this->_hasComputedMapping = false;
176}
177
178template <typename SOLVER_T, typename... Args>
180{
181 if constexpr (std::tuple_size_v<std::tuple<Args...>>> 0) {
182 auto param = std::get<0>(optionalArgs);
183 std::string exec = param.executor;
184 if (param.solver == "qr-solver") {
185 return "global-direct RBF (" + exec + ")";
186 } else {
187 return "global-iterative RBF (" + exec + ")";
188 }
189 } else {
190 return "global-direct RBF (cpu-executor)";
191 }
192}
193
194template <typename SOLVER_T, typename... Args>
196{
198 precice::profiling::Event e("map.rbf.mapData.From" + this->input()->getName() + "To" + this->output()->getName(), profiling::Synchronize);
200
201 PRECICE_DEBUG("Map conservative using {}", getName());
202
203 // Gather input data
205
206 const auto &localInData = inData.values;
207
208 int localOutputSize = 0;
209 for (const auto &vertex : this->output()->vertices()) {
210 if (vertex.isOwner()) {
211 ++localOutputSize;
212 }
213 }
214
215 localOutputSize *= inData.dataDims;
216
217 utils::IntraComm::getCommunication()->sendRange(localInData, 0);
218 utils::IntraComm::getCommunication()->send(localOutputSize, 0);
219
220 } else { // Parallel Primary rank or Serial case
221
222 std::vector<double> globalInValues;
223 std::vector<double> outputValueSizes;
224 {
225 const auto &localInData = inData.values;
226 globalInValues.insert(globalInValues.begin(), localInData.data(), localInData.data() + localInData.size());
227
228 int localOutputSize = 0;
229 for (const auto &vertex : this->output()->vertices()) {
230 if (vertex.isOwner()) {
231 ++localOutputSize;
232 }
233 }
234
235 localOutputSize *= inData.dataDims;
236
237 outputValueSizes.push_back(localOutputSize);
238 }
239
240 {
241 int secondaryOutputValueSize;
243 std::vector<double> secondaryBuffer = utils::IntraComm::getCommunication()->receiveRange(rank, AsVectorTag<double>{});
244 globalInValues.insert(globalInValues.end(), secondaryBuffer.begin(), secondaryBuffer.end());
245
246 utils::IntraComm::getCommunication()->receive(secondaryOutputValueSize, rank);
247 outputValueSizes.push_back(secondaryOutputValueSize);
248 }
249 }
250
251 const int valueDim = inData.dataDims;
252
253 // Construct Eigen vectors
254 Eigen::Map<Eigen::VectorXd> inputValues(globalInValues.data(), globalInValues.size());
255 Eigen::VectorXd outputValues((this->output()->getGlobalNumberOfVertices()) * valueDim);
256
257 Eigen::VectorXd in; // rows == outputSize
258 in.resize(_rbfSolver->getOutputSize()); // rows == outputSize
259
260 outputValues.setZero();
261
262 for (int dim = 0; dim < valueDim; dim++) {
263 for (int i = 0; i < in.size(); i++) { // Fill input data values
264 in[i] = inputValues(i * valueDim + dim);
265 }
266
267 Eigen::VectorXd out;
268 out = _rbfSolver->solveConservative(in, _polynomial);
269
270 // Copy mapped data to output data values
271 for (int i = 0; i < this->output()->getGlobalNumberOfVertices(); i++) {
272 outputValues[i * valueDim + dim] = out[i];
273 }
274 }
275
276 // Data scattering to secondary ranks
278
279 // Filter data
280 int outputCounter = 0;
281 for (int i = 0; i < static_cast<int>(this->output()->nVertices()); ++i) {
282 if (this->output()->vertex(i).isOwner()) {
283 for (int dim = 0; dim < valueDim; ++dim) {
284 outData[i * valueDim + dim] = outputValues(outputCounter);
285 ++outputCounter;
286 }
287 }
288 }
289
290 // Data scattering to secondary ranks
291 int beginPoint = outputValueSizes.at(0);
293 precice::span<const double> toSend{outputValues.data() + beginPoint, static_cast<size_t>(outputValueSizes.at(rank))};
294 utils::IntraComm::getCommunication()->sendRange(toSend, rank);
295 beginPoint += outputValueSizes.at(rank);
296 }
297 } else { // Serial
298 outData = outputValues;
299 }
300 }
302 std::vector<double> receivedValues = utils::IntraComm::getCommunication()->receiveRange(0, AsVectorTag<double>{});
303
304 const int valueDim = inData.dataDims;
305
306 int outputCounter = 0;
307 for (int i = 0; i < static_cast<int>(this->output()->nVertices()); ++i) {
308 if (this->output()->vertex(i).isOwner()) {
309 for (int dim = 0; dim < valueDim; ++dim) {
310 outData[i * valueDim + dim] = receivedValues.at(outputCounter);
311 ++outputCounter;
312 }
313 }
314 }
315 }
316}
317
318template <typename SOLVER_T, typename... Args>
320{
322 precice::profiling::Event e("map.rbf.mapData.From" + this->input()->getName() + "To" + this->output()->getName(), profiling::Synchronize);
324
325 PRECICE_DEBUG("Map {} using {}", (this->hasConstraint(Mapping::CONSISTENT) ? "consistent" : "scaled-consistent"), getName());
326
327 // Gather input data
329 // Input data is filtered
330 auto localInDataFiltered = this->input()->getOwnedVertexData(inData.values);
331 int localOutputSize = outData.size();
332
333 // Send data and output size
334 utils::IntraComm::getCommunication()->sendRange(localInDataFiltered, 0);
335 utils::IntraComm::getCommunication()->send(localOutputSize, 0);
336
337 } else { // Primary rank or Serial case
338
339 const int valueDim = inData.dataDims;
340
341 std::vector<double> globalInValues(static_cast<std::size_t>(this->input()->getGlobalNumberOfVertices()) * valueDim, 0.0);
342 std::vector<int> outValuesSize;
343
344 if (utils::IntraComm::isPrimary()) { // Parallel case
345
346 // Filter input data
347 const auto &localInData = this->input()->getOwnedVertexData(inData.values);
348 std::copy(localInData.data(), localInData.data() + localInData.size(), globalInValues.begin());
349 outValuesSize.push_back(outData.size());
350
351 int inputSizeCounter = localInData.size();
352 int secondaryOutDataSize{0};
353
355 std::vector<double> secondaryBuffer = utils::IntraComm::getCommunication()->receiveRange(rank, AsVectorTag<double>{});
356 std::copy(secondaryBuffer.begin(), secondaryBuffer.end(), globalInValues.begin() + inputSizeCounter);
357 inputSizeCounter += secondaryBuffer.size();
358
359 utils::IntraComm::getCommunication()->receive(secondaryOutDataSize, rank);
360 outValuesSize.push_back(secondaryOutDataSize);
361 }
362
363 } else { // Serial case
364 const auto &localInData = inData.values;
365 std::copy(localInData.data(), localInData.data() + localInData.size(), globalInValues.begin());
366 outValuesSize.push_back(outData.size());
367 }
368
369 Eigen::VectorXd in;
370
371 in.resize(_rbfSolver->getInputSize()); // rows == n
372 in.setZero();
373
374 // Construct Eigen vectors
375 Eigen::Map<Eigen::VectorXd> inputValues(globalInValues.data(), globalInValues.size());
376
377 Eigen::VectorXd outputValues;
378 outputValues.resize((_rbfSolver->getOutputSize()) * valueDim);
379
380 Eigen::VectorXd out;
381 outputValues.setZero();
382
383 // For every data dimension, perform mapping
384 for (int dim = 0; dim < valueDim; dim++) {
385 // Fill input from input data values (last polyparams entries remain zero)
386 for (int i = 0; i < this->input()->getGlobalNumberOfVertices(); i++) {
387 in[i] = inputValues[i * valueDim + dim];
388 }
389
390 out = _rbfSolver->solveConsistent(in, _polynomial);
391
392 // Copy mapped data to output data values
393 for (int i = 0; i < out.size(); i++) {
394 outputValues[i * valueDim + dim] = out[i];
395 }
396 }
397
398 outData = Eigen::Map<Eigen::VectorXd>(outputValues.data(), outValuesSize.at(0));
399
400 // Data scattering to secondary ranks
401 int beginPoint = outValuesSize.at(0);
402
405 precice::span<const double> toSend{outputValues.data() + beginPoint, static_cast<size_t>(outValuesSize.at(rank))};
406 utils::IntraComm::getCommunication()->sendRange(toSend, rank);
407 beginPoint += outValuesSize.at(rank);
408 }
409 }
410 }
412 std::vector<double> receivedValues = utils::IntraComm::getCommunication()->receiveRange(0, AsVectorTag<double>{});
413 outData = Eigen::Map<Eigen::VectorXd>(receivedValues.data(), receivedValues.size());
414 }
415}
416} // namespace mapping
417} // namespace precice
std::ostream & out
#define PRECICE_DEBUG(...)
Definition LogMacros.hpp:64
#define PRECICE_TRACE(...)
Definition LogMacros.hpp:95
#define PRECICE_CHECK(check,...)
Definition LogMacros.hpp:35
#define PRECICE_ASSERT(...)
Definition assertion.hpp:87
T at(T... args)
T begin(T... args)
Abstract base class for mapping of data from one mesh to another.
Definition Mapping.hpp:15
Constraint
Specifies additional constraints for a mapping.
Definition Mapping.hpp:29
InitialGuessRequirement
Specifies whether the mapping requires an initial guess.
Definition Mapping.hpp:63
Mapping with radial basis functions.
typename SOLVER_T::BASIS_FUNCTION_T RADIAL_BASIS_FUNCTION_T
Polynomial _polynomial
Treatment of the polynomial.
std::tuple< Args... > optionalArgs
Optional constructor arguments for the solver class.
void computeMapping() final override
Computes the mapping coefficients from the in- and output mesh.
void mapConsistent(const time::Sample &inData, Eigen::VectorXd &outData) final override
Maps data using a consistent constraint.
void clear() final override
Removes a computed mapping.
RadialBasisFctMapping(Mapping::Constraint constraint, int dimensions, RADIAL_BASIS_FUNCTION_T function, std::array< bool, 3 > deadAxis, Polynomial polynomial, Args... args)
Constructor.
std::string getName() const final override
name of the rbf mapping
void mapConservative(const time::Sample &inData, Eigen::VectorXd &outData) final override
Maps data using a conservative constraint.
Container and creator for meshes.
Definition Mesh.hpp:39
void addMesh(Mesh &deltaMesh)
Definition Mesh.cpp:309
static constexpr MeshID MESH_ID_UNDEFINED
Use if the id of the mesh is not necessary.
Definition Mesh.hpp:58
std::size_t nVertices() const
Returns the number of vertices.
Definition Mesh.cpp:63
Vertex of a mesh.
Definition Vertex.hpp:16
bool isOwner() const
Definition Vertex.cpp:22
A C++ 11 implementation of the non-owning C++20 std::span type.
Definition span.hpp:284
constexpr pointer data() const noexcept
Definition span.hpp:500
static bool isPrimary()
True if this process is running the primary rank.
Definition IntraComm.cpp:52
static auto allSecondaryRanks()
Returns an iterable range over salve ranks [1, _size)
Definition IntraComm.hpp:37
static bool isSecondary()
True if this process is running a secondary rank.
Definition IntraComm.cpp:57
static com::PtrCommunication & getCommunication()
Intra-participant communication.
Definition IntraComm.hpp:31
T copy(T... args)
T data(T... args)
T end(T... args)
T insert(T... args)
void sendMesh(Communication &communication, int rankReceiver, const mesh::Mesh &mesh)
Definition Extra.cpp:8
void receiveMesh(Communication &communication, int rankSender, mesh::Mesh &mesh)
Definition Extra.cpp:13
Polynomial
How to handle the polynomial?
void filterMesh(Mesh &destination, const Mesh &source, UnaryPredicate p)
Definition Filter.hpp:17
static constexpr SynchronizeTag Synchronize
Convenience instance of the SynchronizeTag.
Definition Event.hpp:21
Main namespace of the precice library.
int Rank
Definition Types.hpp:37
STL namespace.
T push_back(T... args)
T size(T... args)
int dataDims
The dimensionality of the data.
Definition Sample.hpp:45
Eigen::VectorXd values
Definition Sample.hpp:49
T tuple_size_v