preCICE v3.1.1
Loading...
Searching...
No Matches
PetRadialBasisFctMapping.hpp
Go to the documentation of this file.
1#pragma once
2#ifndef PRECICE_NO_PETSC
3
4#include <iterator>
5#include <map>
6#include <numeric>
7#include <vector>
8
10#include "logging/LogMacros.hpp"
13#include "math/math.hpp"
14#include "precice/impl/versions.hpp"
15#include "profiling/Event.hpp"
16#include "utils/IntraComm.hpp"
17#include "utils/Petsc.hpp"
18#include "utils/assertion.hpp"
19
21
22namespace precice {
23namespace mapping {
24
25namespace {
26// VecChop was deprecated in PETSc 3.20 and is to be replaced by VecFilter
27inline PetscErrorCode PRECICE_VecFilter(Vec v, PetscReal tol)
28{
29#if ((PETSC_MAJOR > 3) || (PETSC_MAJOR == 3 && PETSC_MINOR >= 20))
30 return VecFilter(v, tol);
31#else
32 return VecChop(v, tol);
33#endif
34}
35} // namespace
36
37namespace tests {
38class PetRadialBasisFctMappingTest; // Forward declaration to friend the class
39}
40
51template <typename RADIAL_BASIS_FUNCTION_T>
52class PetRadialBasisFctMapping : public RadialBasisFctBaseMapping<RADIAL_BASIS_FUNCTION_T> {
53public:
67 Mapping::Constraint constraint,
68 int dimensions,
69 const RADIAL_BASIS_FUNCTION_T &function,
70 std::array<bool, 3> deadAxis,
71 double solverRtol = 1e-9,
72 Polynomial polynomial = Polynomial::SEPARATE);
73
76
78 void computeMapping() final override;
79
81 void clear() final override;
82
84 std::string getName() const final override;
85
86private:
88 void mapConservative(const time::Sample &inData, Eigen::VectorXd &outData) final override;
89
91 void mapConsistent(const time::Sample &inData, Eigen::VectorXd &outData) final override;
92
95
96 mutable logging::Logger _log{"mapping::PetRadialBasisFctMapping"};
97
100
103
106
109
112
115
118
121
122 const double _solverRtol;
123
126
128 bool useRescaling = true;
129
132
135
138
140 void printMappingInfo(int dim) const;
141
144
147
149
156 void loadInitialGuessForDim(int dimension, int allDimensions, petsc::Vector &destination);
157
164 void storeInitialGuessForDim(int dimension, int allDimensions, petsc::Vector &source);
165};
166
167// --------------------------------------------------- HEADER IMPLEMENTATIONS
168
169template <typename RADIAL_BASIS_FUNCTION_T>
171 Mapping::Constraint constraint,
172 int dimensions,
173 const RADIAL_BASIS_FUNCTION_T &function,
174 std::array<bool, 3> deadAxis,
175 double solverRtol,
176 Polynomial polynomial)
177 : RadialBasisFctBaseMapping<RADIAL_BASIS_FUNCTION_T>(constraint, dimensions, function, deadAxis, Mapping::InitialGuessRequirement::Required),
178 _matrixC("C"),
179 _matrixQ("Q"),
180 _matrixA("A"),
181 _matrixV("V"),
182 _solver("Coefficient Solver"),
183 _QRsolver("QR Solver"),
184 _AOmapping(nullptr),
185 _solverRtol(solverRtol),
186 _polynomial(polynomial),
187 _commState(utils::Parallel::current())
188{
191 localPolyparams = (_commState->rank() > 0) ? 0 : polyparams;
192}
193
194template <typename RADIAL_BASIS_FUNCTION_T>
199
200template <typename RADIAL_BASIS_FUNCTION_T>
202{
204 precice::profiling::Event e("map.pet.computeMapping.From" + this->input()->getName() + "To" + this->output()->getName(), profiling::Synchronize);
205 precice::profiling::Event ePreCompute("map.pet.preComputeMapping.From" + this->input()->getName() + "To" + this->output()->getName());
206
207 clear();
208
209 PRECICE_DEBUG_IF(_polynomial == Polynomial::ON, "Using integrated polynomial.");
210 PRECICE_DEBUG_IF(_polynomial == Polynomial::OFF, "Using no polynomial.");
211 PRECICE_DEBUG_IF(_polynomial == Polynomial::SEPARATE, "Using separated polynomial.");
212
213 PRECICE_ASSERT(this->input()->getDimensions() == this->output()->getDimensions(),
214 this->input()->getDimensions(), this->output()->getDimensions());
215 int const dimensions = this->input()->getDimensions();
216 mesh::PtrMesh inMesh;
217 mesh::PtrMesh outMesh;
218 if (this->hasConstraint(Mapping::CONSERVATIVE)) {
219 inMesh = this->output();
220 outMesh = this->input();
221 } else {
222 inMesh = this->input();
223 outMesh = this->output();
224 }
225
226 // Indizes that are used to build the Petsc AO mapping
227 std::vector<PetscInt> myIndizes;
228
229 // Indizes for Q^T, holding the polynomial
230 if (_commState->rank() <= 0) // Rank 0 or not in IntraComm mode
231 for (size_t i = 0; i < polyparams; i++)
232 myIndizes.push_back(i); // polyparams reside in the first rows (which are always on rank 0)
233
234 // Indizes for the vertices with polyparams offset
235 for (const mesh::Vertex &v : inMesh->vertices())
236 if (v.isOwner())
237 myIndizes.push_back(v.getGlobalIndex() + polyparams);
238
239 auto n = myIndizes.size(); // polyparams, if on rank 0, are included here
240
241 auto outputSize = outMesh->nVertices();
242
243 PetscErrorCode ierr = 0;
244 PRECICE_DEBUG("inMesh->nVertices() = {}", inMesh->nVertices());
245 PRECICE_DEBUG("outMesh->nVertices() = {}", outMesh->nVertices());
246 ePreCompute.stop();
247
248 precice::profiling::Event eCreateMatrices("map.pet.createMatrices.From" + this->input()->getName() + "To" + this->output()->getName(), profiling::Synchronize);
249
250 // Matrix C: Symmetric, sparse matrix with n x n local size.
251 _matrixC.init(n, n, PETSC_DETERMINE, PETSC_DETERMINE, MATSBAIJ);
252 PRECICE_DEBUG("Set matrix C to local size {0} x {0}", n);
253 ierr = MatSetOption(_matrixC, MAT_SYMMETRIC, PETSC_TRUE);
254 CHKERRV(ierr);
255 ierr = MatSetOption(_matrixC, MAT_SYMMETRY_ETERNAL, PETSC_TRUE);
256 CHKERRV(ierr);
257
258 // Matrix Q: Dense, holds the input mesh for the polynomial if set to SEPARATE. Zero size otherwise
259 _matrixQ.init(n, PETSC_DETERMINE, PETSC_DETERMINE, sepPolyparams, MATDENSE);
260 PRECICE_DEBUG("Set matrix Q to local size {} x {}", n, sepPolyparams);
261
262 // Matrix V: Dense, holds the output mesh for polynomial if set to SEPARATE. Zero size otherwise
263 _matrixV.init(outputSize, PETSC_DETERMINE, PETSC_DETERMINE, sepPolyparams, MATDENSE);
264 PRECICE_DEBUG("Set matrix V to local size {} x {}", outputSize, sepPolyparams);
265
266 // Matrix A: Sparse matrix with outputSize x n local size.
267 _matrixA.init(outputSize, n, PETSC_DETERMINE, PETSC_DETERMINE, MATAIJ);
268 PRECICE_DEBUG("Set matrix A to local size {} x {}", outputSize, n);
269
270 eCreateMatrices.stop();
271 precice::profiling::Event eAO("map.pet.AO.From" + this->input()->getName() + "To" + this->output()->getName(), profiling::Synchronize);
272
273 auto const ownerRangeABegin = _matrixA.ownerRange().first;
274 auto const ownerRangeAEnd = _matrixA.ownerRange().second;
275
276 // A mapping from globalIndex -> local col/row
277 ierr = AOCreateMapping(_commState->comm, myIndizes.size(), myIndizes.data(), nullptr, &_AOmapping);
278 CHKERRV(ierr);
279
280 eAO.stop();
281
282 Eigen::VectorXd distance(dimensions);
283
284 // We do preallocating of the matrices C and A. That means we traverse the input data once, just
285 // to know where we have entries in the sparse matrix. This information petsc can use to
286 // preallocate the matrix. In the second phase we actually fill the matrix.
287
288 // Stores col -> value for each row;
289 VertexData vertexData = bgPreallocationMatrixC(inMesh);
290
291 // -- BEGIN FILL LOOP FOR MATRIX C --
292 PRECICE_DEBUG("Begin filling matrix C");
293 precice::profiling::Event eFillC("map.pet.fillC.From" + this->input()->getName() + "To" + this->output()->getName(), profiling::Synchronize);
294
295 // We collect entries for each row and set them blockwise using MatSetValues.
296 PetscInt const idxSize = std::max(_matrixC.getSize().second, _matrixQ.getSize().second);
297 std::vector<PetscInt> colIdx(idxSize); // holds the columns indices of the entries
298 std::vector<PetscScalar> rowVals(idxSize); // holds the values of the entries
299
300 int preallocRow = 0;
301 PetscInt row = _matrixC.ownerRange().first + localPolyparams;
302 for (const mesh::Vertex &inVertex : inMesh->vertices()) {
303 if (not inVertex.isOwner())
304 continue;
305
306 PetscInt colNum = 0; // holds the number of non-zero columns in current row
307
308 // -- SETS THE POLYNOMIAL PART OF THE MATRIX --
309 if (_polynomial == Polynomial::ON or _polynomial == Polynomial::SEPARATE) {
310 colIdx[colNum] = colNum;
311 rowVals[colNum++] = 1;
312
313 for (int dim = 0; dim < dimensions; dim++) {
314 if (not this->_deadAxis[dim]) {
315 colIdx[colNum] = colNum;
316 rowVals[colNum++] = inVertex.coord(dim);
317 }
318 }
319
320 // cols are always the first ones for the polynomial, no need to translate
321 if (_polynomial == Polynomial::ON) {
322 ierr = MatSetValues(_matrixC, colNum, colIdx.data(), 1, &row, rowVals.data(), INSERT_VALUES);
323 CHKERRV(ierr);
324 } else if (_polynomial == Polynomial::SEPARATE) {
325 ierr = MatSetValues(_matrixQ, 1, &row, colNum, colIdx.data(), rowVals.data(), INSERT_VALUES);
326 CHKERRV(ierr);
327 }
328 colNum = 0;
329 }
330
331 // -- SETS THE COEFFICIENTS --
332 {
333 auto const &rowVertices = vertexData[preallocRow];
334 for (const auto &vertex : rowVertices) {
335 rowVals[colNum] = this->_basisFunction.evaluate(vertex.second);
336 colIdx[colNum++] = vertex.first;
337 }
338 ++preallocRow;
339 }
340 ierr = AOApplicationToPetsc(_AOmapping, colNum, colIdx.data());
341 CHKERRV(ierr);
342 ierr = MatSetValues(_matrixC, 1, &row, colNum, colIdx.data(), rowVals.data(), INSERT_VALUES);
343 CHKERRV(ierr);
344 ++row;
345 }
346 PRECICE_DEBUG("Finished filling Matrix C");
347 eFillC.stop();
348 // -- END FILL LOOP FOR MATRIX C --
349
350 // PETSc requires that all diagonal entries are set, even if set to zero.
351 _matrixC.assemble(MAT_FLUSH_ASSEMBLY);
352 auto zeros = petsc::Vector::allocate(_matrixC);
353 VecZeroEntries(zeros);
354 zeros.assemble();
355 //MatDiagonalSet(_matrixC, zeros, INSERT_VALUES);
356 MatDiagonalSet(_matrixC, zeros, ADD_VALUES);
357
358 // Begin assembly here, all assembly is ended at the end of this function.
359 ierr = MatAssemblyBegin(_matrixC, MAT_FINAL_ASSEMBLY);
360 CHKERRV(ierr);
361 ierr = MatAssemblyBegin(_matrixQ, MAT_FINAL_ASSEMBLY);
362 CHKERRV(ierr);
363
364 vertexData = bgPreallocationMatrixA(inMesh, outMesh);
365
366 // holds the columns indices of the entries
367 colIdx.resize(std::max(_matrixA.getSize().second, _matrixV.getSize().second));
368 // holds the values of the entries
369 rowVals.resize(std::max(_matrixA.getSize().second, _matrixV.getSize().second));
370
371 // -- BEGIN FILL LOOP FOR MATRIX A --
372 PRECICE_DEBUG("Begin filling matrix A.");
373 precice::profiling::Event eFillA("map.pet.fillA.From" + this->input()->getName() + "To" + this->output()->getName(), profiling::Synchronize);
374
375 for (PetscInt row = ownerRangeABegin; row < ownerRangeAEnd; ++row) {
376 mesh::Vertex const &oVertex = outMesh->vertex(row - _matrixA.ownerRange().first);
377
378 // -- SET THE POLYNOMIAL PART OF THE MATRIX --
379 if (_polynomial == Polynomial::ON or _polynomial == Polynomial::SEPARATE) {
380 petsc::Matrix *m = _polynomial == Polynomial::ON ? &_matrixA : &_matrixV;
381 PetscInt colNum = 0;
382
383 colIdx[colNum] = colNum;
384 rowVals[colNum++] = 1;
385
386 for (int dim = 0; dim < dimensions; dim++) {
387 if (not this->_deadAxis[dim]) {
388 colIdx[colNum] = colNum;
389 rowVals[colNum++] = oVertex.coord(dim);
390 }
391 }
392 ierr = MatSetValues(*m, 1, &row, colNum, colIdx.data(), rowVals.data(), INSERT_VALUES);
393 CHKERRV(ierr);
394 }
395
396 // -- SETS THE COEFFICIENTS --
397 PetscInt colNum = 0;
398
399 {
400 auto const &rowVertices = vertexData[row - ownerRangeABegin];
401 for (const auto &vertex : rowVertices) {
402 rowVals[colNum] = this->_basisFunction.evaluate(vertex.second);
403 colIdx[colNum++] = vertex.first;
404 }
405 }
406 ierr = AOApplicationToPetsc(_AOmapping, colNum, colIdx.data());
407 CHKERRV(ierr);
408 ierr = MatSetValues(_matrixA, 1, &row, colNum, colIdx.data(), rowVals.data(), INSERT_VALUES);
409 CHKERRV(ierr);
410 }
411 PRECICE_DEBUG("Finished filling Matrix A");
412 eFillA.stop();
413 // -- END FILL LOOP FOR MATRIX A --
414
415 precice::profiling::Event ePostFill("map.pet.postFill.From" + this->input()->getName() + "To" + this->output()->getName(), profiling::Synchronize);
416
417 ierr = MatAssemblyBegin(_matrixA, MAT_FINAL_ASSEMBLY);
418 CHKERRV(ierr);
419
420 ierr = MatAssemblyEnd(_matrixC, MAT_FINAL_ASSEMBLY);
421 CHKERRV(ierr);
422 ierr = MatAssemblyEnd(_matrixQ, MAT_FINAL_ASSEMBLY);
423 CHKERRV(ierr);
424 ierr = MatAssemblyEnd(_matrixA, MAT_FINAL_ASSEMBLY);
425 CHKERRV(ierr);
426 _matrixV.assemble();
427
428 ePostFill.stop();
429
430 precice::profiling::Event eSolverInit("map.pet.solverInit.From" + this->input()->getName() + "To" + this->output()->getName(), profiling::Synchronize);
431
432 // -- CONFIGURE SOLVER FOR POLYNOMIAL --
433 if (_polynomial == Polynomial::SEPARATE) {
434 PC pc;
435 KSPGetPC(_QRsolver, &pc);
436 PCSetType(pc, PCNONE);
437 KSPSetType(_QRsolver, KSPLSQR);
438 KSPSetOperators(_QRsolver, _matrixQ, _matrixQ);
439 }
440
441 // -- CONFIGURE SOLVER FOR SYSTEM MATRIX --
442 KSPSetOperators(_solver, _matrixC, _matrixC);
443 CHKERRV(ierr);
444 // The fourth argument defines the divergence tolerance, i.e. the tolerance, for which the linear system is
445 // considered as diverged and the computation is aborted. The PETSC_DEFAULT value is 1e9. However, the
446 // divergence residual is scaled using the RHS b of the system. In case the RHS is (still nonzero) very small
447 // and we have a (bad) nonzero initial guess x as defined below, the divergence tolerance might be exceeded
448 // in the first iteration, although the system could be solved in the usual way. This behavior is essentially
449 // a glitch in the divergence check. Therefore, we select a very high value (1e30) in order to disable the
450 // divergence check. In practice, the check is very rarely needed with Krylov methods. According to the PETSc
451 // people, the rare use cases aim for a bad preconditioner, which is not even used in our configuration.
452 // Hence, we can disable the divergence check without concerns.
453 KSPSetTolerances(_solver, _solverRtol, PETSC_DEFAULT, 1e30, PETSC_DEFAULT);
454 KSPSetInitialGuessNonzero(_solver, PETSC_TRUE);
455 CHKERRV(ierr); // Reuse the results from the last iteration, held in the out vector.
456 KSPSetOptionsPrefix(_solver, "solverC_"); // s.t. options for only this solver can be set on the command line
457 KSPSetFromOptions(_solver);
458
459 eSolverInit.stop();
460
461 // if (totalNNZ > static_cast<size_t>(20*n)) {
462 // PRECICE_DEBUG("Using Cholesky decomposition as direct solver for dense matrix.");
463 // PC prec;
464 // KSPSetType(_solver, KSPPREONLY);
465 // KSPGetPC(_solver, &prec);
466 // PCSetType(prec, PCCHOLESKY);
467 // PCFactorSetShiftType(prec, MAT_SHIFT_NONZERO);
468 // }
469
470 // -- COMPUTE RESCALING COEFFICIENTS USING THE SYSTEM MATRIX C SOLVER --
471 if (useRescaling and (_polynomial == Polynomial::SEPARATE)) {
472 precice::profiling::Event eRescaling("map.pet.computeRescaling.From" + this->input()->getName() + "To" + this->output()->getName(), profiling::Synchronize);
473 auto rhs = petsc::Vector::allocate(_matrixC);
474 auto rescalingCoeffs = petsc::Vector::allocate(_matrixC);
475 VecSet(rhs, 1);
476 rhs.assemble();
477 if (_solver.solve(rhs, rescalingCoeffs) == petsc::KSPSolver::SolverResult::Converged) {
478 PRECICE_INFO("Using rescaling. {}", _solver.summaryFor(rhs));
479 } else {
480 PRECICE_WARN("Deactivating rescaling! {}", _solver.summaryFor(rhs));
481 useRescaling = false;
482 }
483
484 eRescaling.addData("Iterations", _solver.getIterationNumber());
485 ierr = MatCreateVecs(_matrixA, nullptr, &oneInterpolant.vector);
486 CHKERRV(ierr);
487 ierr = MatMult(_matrixA, rescalingCoeffs, oneInterpolant);
488 CHKERRV(ierr); // get the output of g(x) = 1
489 // set values close to zero to exactly 0.0, s.t. PointwiseDevide doesn't do division on these entries
490 ierr = PRECICE_VecFilter(oneInterpolant, 1e-6);
491 CHKERRV(ierr);
492 }
493
494 this->_hasComputedMapping = true;
495
496 PRECICE_DEBUG("Number of mallocs for matrix C = {}", _matrixC.getInfo(MAT_LOCAL).mallocs);
497 PRECICE_DEBUG("Non-zeros allocated / used / unused for matrix C = {} / {} / {}",
498 _matrixC.getInfo(MAT_LOCAL).nz_allocated, _matrixC.getInfo(MAT_LOCAL).nz_used, _matrixC.getInfo(MAT_LOCAL).nz_unneeded);
499 PRECICE_DEBUG("Number of mallocs for matrix A = {}", _matrixA.getInfo(MAT_LOCAL).mallocs);
500 PRECICE_DEBUG("Non-zeros allocated / used / unused for matrix A = {} / {} / {}",
501 _matrixA.getInfo(MAT_LOCAL).nz_allocated, _matrixA.getInfo(MAT_LOCAL).nz_used, _matrixA.getInfo(MAT_LOCAL).nz_unneeded);
502}
503
504template <typename RADIAL_BASIS_FUNCTION_T>
506{
507 _matrixC.reset();
508 _matrixA.reset();
509 _matrixQ.reset();
510 _matrixV.reset();
511
512 _solver.reset();
513 _QRsolver.reset();
514
515 petsc::destroy(&_AOmapping);
516
517 this->_hasComputedMapping = false;
518}
519
520template <typename RADIAL_BASIS_FUNCTION_T>
522{
523 return "global-iterative RBF";
524}
525
526template <typename RADIAL_BASIS_FUNCTION_T>
528{
529 // Only skip collectively over all MPI ranks as we use collective Vector ops
530 if (destination.getSize() == 0) {
531 return;
532 }
533 auto sizePerDim = destination.getLocalSize();
534 auto totalSize = sizePerDim * allDimensions;
535
536 if (!this->hasInitialGuess() && (totalSize > 0)) {
537 // We don't need to modify the petsc vector here
538 this->initialGuess() = Eigen::VectorXd::Zero(totalSize);
539 }
540
541 PRECICE_ASSERT(this->initialGuess().size() == totalSize, this->initialGuess().size(), totalSize);
542 auto offset = dimension * sizePerDim;
543 auto begin = std::next(this->initialGuess().data(), offset);
544 destination.copyFrom({begin, static_cast<std::size_t>(sizePerDim)});
545}
546
547template <typename RADIAL_BASIS_FUNCTION_T>
549{
550 // Only skip collectively over all MPI ranks as we use collective Vector ops
551 if (source.getSize() == 0) {
552 return;
553 }
554 auto sizePerDim = source.getLocalSize();
555
556 PRECICE_ASSERT(this->hasInitialGuess() || (sizePerDim == 0), "Call loadInitialGuessForDim first");
557 PRECICE_ASSERT(this->initialGuess().size() == sizePerDim * allDimensions, this->initialGuess().size(), sizePerDim * allDimensions);
558 auto offset = dimension * sizePerDim;
559 auto begin = std::next(this->initialGuess().data(), offset);
560 source.copyTo({begin, static_cast<std::size_t>(sizePerDim)});
561}
562
563template <typename RADIAL_BASIS_FUNCTION_T>
565{
567 precice::profiling::Event e("map.pet.mapData.From" + this->input()->getName() + "To" + this->output()->getName(), profiling::Synchronize);
568
569 PetscErrorCode ierr = 0;
570 auto const & inValues = inData.values;
571 auto & outValues = outData;
572
573 int const valueDim = inData.dataDims;
574 PRECICE_ASSERT(this->hasConstraint(Mapping::CONSISTENT) || this->isScaledConsistent());
575
576 auto out = petsc::Vector::allocate(_matrixA, "out");
577 auto in = petsc::Vector::allocate(_matrixC, "in");
578 auto a = petsc::Vector::allocate(_matrixQ, "a", petsc::Vector::RIGHT); // holds the solution of the LS polynomial
579
580 PetscScalar const *vecArray;
581
582 // For every data dimension, perform mapping
583 for (int dim = 0; dim < valueDim; dim++) {
584 printMappingInfo(dim);
585
586 // Fill input from input data values
588 inVals.reserve(this->input()->nVertices());
590 inIdx.reserve(this->input()->nVertices());
591 for (size_t i = 0; i < this->input()->nVertices(); ++i) {
592 if (not this->input()->vertex(i).isOwner())
593 continue;
594 inVals.emplace_back(inValues[i * valueDim + dim]);
595 inIdx.emplace_back(inIdx.size() + in.ownerRange().first + localPolyparams);
596 }
597 ierr = VecSetValues(in, inIdx.size(), inIdx.data(), inVals.data(), INSERT_VALUES);
598 CHKERRV(ierr);
599 in.assemble();
600
601 if (_polynomial == Polynomial::SEPARATE) {
602 switch (_QRsolver.solve(in, a)) {
603 case (petsc::KSPSolver::SolverResult::Converged):
604 PRECICE_DEBUG("The polynomial QR system of the RBF mapping from mesh {} to mesh {} converged. {}",
605 this->input()->getName(), this->output()->getName(), _QRsolver.summaryFor(in));
606 break;
607 case (petsc::KSPSolver::SolverResult::Stopped):
608 PRECICE_WARN("The polynomial QR system of the RBF mapping from mesh {} to mesh {} has not converged. "
609 "This means most probably that the mapping problem is not well-posed or your relative tolerance is too conservative. "
610 "Please check if your coupling meshes are correct. "
611 "Maybe you need to fix axis-aligned mapping setups by marking perpendicular axes as dead? {}",
612 this->input()->getName(), this->output()->getName(), _QRsolver.summaryFor(in));
613 break;
614 case (petsc::KSPSolver::SolverResult::Diverged):
615 KSPView(_QRsolver, PETSC_VIEWER_STDOUT_WORLD);
616 PRECICE_ERROR("The polynomial QR system of the RBF mapping from mesh {} to mesh {} has diverged. "
617 "This means most probably that the mapping problem is not well-posed. "
618 "Please check if your coupling meshes are correct. "
619 "Maybe you need to fix axis-aligned mapping setups by marking perpendicular axes as dead? {}",
620 this->input()->getName(), this->output()->getName(), _QRsolver.summaryFor(in));
621 break;
622 }
623
624 VecScale(a, -1);
625 // in = Q * a + in
626 MatMultAdd(_matrixQ, a, in, in); // Subtract the polynomial from the input values
627 }
628
629 petsc::Vector p = petsc::Vector::allocate(_matrixC, "p");
630 loadInitialGuessForDim(dim, valueDim, p);
631
632 profiling::Event eSolve("map.pet.solveConsistent.From" + this->input()->getName() + "To" + this->output()->getName(), profiling::Synchronize);
633 const auto solverResult = _solver.solve(in, p);
634 eSolve.addData("Iterations", _solver.getIterationNumber());
635 eSolve.stop();
636
637 storeInitialGuessForDim(dim, valueDim, p);
638
639 switch (solverResult) {
640 case (petsc::KSPSolver::SolverResult::Converged):
641 PRECICE_DEBUG("The linear system of the RBF mapping from mesh {} to mesh {} converged. {}",
642 this->input()->getName(), this->output()->getName(), _solver.summaryFor(in));
643 break;
644 case (petsc::KSPSolver::SolverResult::Stopped):
645 PRECICE_WARN("The linear system of the RBF mapping from mesh {} to mesh {} has not converged. "
646 "This means most probably that the mapping problem is not well-posed or your relative tolerance is too conservative. "
647 "Please check if your coupling meshes are correct. "
648 "Maybe you need to fix axis-aligned mapping setups by marking perpendicular axes as dead? {}",
649 this->input()->getName(), this->output()->getName(), _solver.summaryFor(in));
650 break;
651 case (petsc::KSPSolver::SolverResult::Diverged):
652 KSPView(_solver, PETSC_VIEWER_STDOUT_WORLD);
653 PRECICE_ERROR("The linear system of the RBF mapping from mesh {} to mesh {} has diverged. "
654 "This means most probably that the mapping problem is not well-posed. "
655 "Please check if your coupling meshes are correct. "
656 "Maybe you need to fix axis-aligned mapping setups by marking perpendicular axes as dead? {}",
657 this->input()->getName(), this->output()->getName(), _solver.summaryFor(in));
658 break;
659 }
660
661 ierr = MatMult(_matrixA, p, out);
662 CHKERRV(ierr);
663
664 if (useRescaling and _polynomial == Polynomial::SEPARATE) {
665 ierr = VecPointwiseDivide(out, out, oneInterpolant);
666 CHKERRV(ierr);
667 }
668
669 if (_polynomial == Polynomial::SEPARATE) {
670 // scale it back to add the polynomial
671 ierr = VecScale(a, -1);
672 // out = V * a + out
673 ierr = MatMultAdd(_matrixV, a, out, out);
674 CHKERRV(ierr);
675 }
676 PRECICE_VecFilter(out, 1e-9);
677
678 // Copy mapped data to output data values
679 ierr = VecGetArrayRead(out, &vecArray);
680 CHKERRV(ierr);
681 int size = out.getLocalSize();
682 for (int i = 0; i < size; i++) {
683 outValues[i * valueDim + dim] = vecArray[i];
684 }
685
686 VecRestoreArrayRead(out, &vecArray);
687 }
688}
689
690template <typename RADIAL_BASIS_FUNCTION_T>
692{
694 precice::profiling::Event e("map.pet.mapData.From" + this->input()->getName() + "To" + this->output()->getName(), profiling::Synchronize);
695
696 PetscErrorCode ierr = 0;
697 auto const & inValues = inData.values;
698 auto & outValues = outData;
699
700 int const valueDim = inData.dataDims;
701 PRECICE_ASSERT(this->hasConstraint(Mapping::CONSERVATIVE));
702
703 auto au = petsc::Vector::allocate(_matrixA, "au", petsc::Vector::RIGHT);
704 auto in = petsc::Vector::allocate(_matrixA, "in");
705 int inRangeStart, inRangeEnd;
706 std::tie(inRangeStart, inRangeEnd) = in.ownerRange();
707 for (int dim = 0; dim < valueDim; dim++) {
708 printMappingInfo(dim);
709
710 // Fill input from input data values
711 for (size_t i = 0; i < this->input()->nVertices(); i++) {
712 auto const globalIndex = this->input()->vertex(i).getGlobalIndex(); // globalIndex is target row
713 if (globalIndex >= inRangeStart and globalIndex < inRangeEnd) // only fill local rows
714 ierr = VecSetValue(in, globalIndex, inValues[i * valueDim + dim], INSERT_VALUES);
715 CHKERRV(ierr);
716 }
717 in.assemble();
718
719 // Gets the petsc::vector for the given combination of outputData, inputData and dimension
720 petsc::Vector out = petsc::Vector::allocate(_matrixC, "out");
721
722 if (_polynomial == Polynomial::SEPARATE) {
723 auto epsilon = petsc::Vector::allocate(_matrixV, "epsilon", petsc::Vector::RIGHT);
724 // epsilon = V^T * in
725 ierr = MatMultTranspose(_matrixV, in, epsilon);
726 CHKERRV(ierr);
727 auto eta = petsc::Vector::allocate(_matrixA, "eta", petsc::Vector::RIGHT);
728 ierr = MatMultTranspose(_matrixA, in, eta);
729 CHKERRV(ierr);
730 auto mu = petsc::Vector::allocate(_matrixC, "mu", petsc::Vector::LEFT);
731 loadInitialGuessForDim(dim, valueDim, mu);
732 _solver.solve(eta, mu);
733 storeInitialGuessForDim(dim, valueDim, mu);
734 VecScale(epsilon, -1);
735 auto tau = petsc::Vector::allocate(_matrixQ, "tau", petsc::Vector::RIGHT);
736 // tau = Q^T * mu + epsilon
737 ierr = MatMultTransposeAdd(_matrixQ, mu, epsilon, tau);
738 CHKERRV(ierr);
739 auto sigma = petsc::Vector::allocate(_matrixQ, "sigma", petsc::Vector::LEFT);
740
741 switch (_QRsolver.solveTranspose(tau, sigma)) {
742 case (petsc::KSPSolver::SolverResult::Converged):
743 PRECICE_DEBUG("The polynomial linear system of the RBF mapping from mesh {} to mesh {} converged. {}",
744 this->input()->getName(), this->output()->getName(), _QRsolver.summaryFor(tau));
745 break;
746 case (petsc::KSPSolver::SolverResult::Stopped):
747 PRECICE_WARN("The polynomial linear system of the RBF mapping from mesh {} to mesh {} has not converged. "
748 "This means most probably that the mapping problem is not well-posed or your relative tolerance is too conservative. "
749 "Please check if your coupling meshes are correct. "
750 "Maybe you need to fix axis-aligned mapping setups by marking perpendicular axes as dead? {}",
751 this->input()->getName(), this->output()->getName(), _QRsolver.summaryFor(tau));
752 break;
753 case (petsc::KSPSolver::SolverResult::Diverged):
754 KSPView(_QRsolver, PETSC_VIEWER_STDOUT_WORLD);
755 PRECICE_ERROR("The polynomial linear system of the RBF mapping from mesh {} to mesh {} "
756 "has diverged. This means most probably that the mapping problem is not well-posed. "
757 "Please check if your coupling meshes are correct. "
758 "Maybe you need to fix axis-aligned mapping setups by marking perpendicular axes as dead? {}",
759 this->input()->getName(), this->output()->getName(), _QRsolver.summaryFor(tau));
760 break;
761 }
762 // out = alpha * sigma + mu.
763 VecWAXPY(out, -1, sigma, mu);
764 } else {
765 ierr = MatMultTranspose(_matrixA, in, au);
766 CHKERRV(ierr);
767
768 loadInitialGuessForDim(dim, valueDim, out);
769
770 profiling::Event eSolve("map.pet.solveConservative.From" + this->input()->getName() + "To" + this->output()->getName(), profiling::Synchronize);
771 const auto solverResult = _solver.solve(au, out);
772 eSolve.addData("Iterations", _solver.getIterationNumber());
773 eSolve.stop();
774
775 storeInitialGuessForDim(dim, valueDim, out);
776
777 switch (solverResult) {
778 case (petsc::KSPSolver::SolverResult::Converged):
779 PRECICE_DEBUG("The linear system of the RBF mapping from mesh {} to mesh {} converged. {}",
780 this->input()->getName(), this->output()->getName(), _solver.summaryFor(au));
781 break;
782 case (petsc::KSPSolver::SolverResult::Stopped):
783 PRECICE_WARN("The linear system of the RBF mapping from mesh {} to mesh {} has not converged. "
784 "This means most probably that the mapping problem is not well-posed or your relative tolerance is too conservative. "
785 "Please check if your coupling meshes are correct. "
786 "Maybe you need to fix axis-aligned mapping setups by marking perpendicular axes as dead? {}",
787 this->input()->getName(), this->output()->getName(), _solver.summaryFor(au));
788 break;
789 case (petsc::KSPSolver::SolverResult::Diverged):
790 KSPView(_solver, PETSC_VIEWER_STDOUT_WORLD);
791 PRECICE_ERROR("The linear system of the RBF mapping from mesh {} to mesh {} "
792 "has diverged. This means most probably that the mapping problem is not well-posed. "
793 "Please check if your coupling meshes are correct. "
794 "Maybe you need to fix axis-aligned mapping setups by marking perpendicular axes as dead? {}",
795 this->input()->getName(), this->output()->getName(), _solver.summaryFor(au));
796 break;
797 }
798 }
799
800 PRECICE_VecFilter(out, 1e-9);
801
802 // Copy mapped data to output data values
803 const PetscScalar *outArray;
804 VecGetArrayRead(out, &outArray);
805
806 int count = 0, ownerCount = 0;
807 for (const mesh::Vertex &vertex : this->output()->vertices()) {
808 if (vertex.isOwner()) {
809 outValues[count * valueDim + dim] = outArray[ownerCount + localPolyparams];
810 ownerCount++;
811 } else {
812 PRECICE_ASSERT(outValues[count * valueDim + dim] == 0.0);
813 }
814 count++;
815 }
816 VecRestoreArrayRead(out, &outArray);
817 }
818}
819
820template <typename RADIAL_BASIS_FUNCTION_T>
822{
823 std::string constraintName;
824 if (this->hasConstraint(Mapping::CONSISTENT)) {
825 constraintName = "consistent";
826 } else if (this->hasConstraint(Mapping::SCALED_CONSISTENT_SURFACE)) {
827 constraintName = "scaled-consistent-surface";
828 } else if (this->hasConstraint(Mapping::SCALED_CONSISTENT_VOLUME)) {
829 constraintName = "scaled-consistent-volume";
830 } else {
831 constraintName = "conservative";
832 }
833
834 const std::string polynomialName = _polynomial == Polynomial::ON ? "on" : _polynomial == Polynomial::OFF ? "off" : "separate";
835
836 PRECICE_INFO("Mapping {} for dimension {} with polynomial set to {}",
837 constraintName, dim, polynomialName);
838}
839
840template <typename RADIAL_BASIS_FUNCTION_T>
843{
844 PRECICE_INFO("Using tree-based preallocation for matrix C");
845 precice::profiling::Event ePreallocC("map.pet.preallocC.From" + this->input()->getName() + "To" + this->output()->getName(), profiling::Synchronize);
846
847 PetscInt n;
848 std::tie(n, std::ignore) = _matrixC.getLocalSize();
849 std::vector<PetscInt> d_nnz(n), o_nnz(n);
850
851 const double supportRadius = this->_basisFunction.getSupportRadius();
852
853 const int dimensions = this->input()->getDimensions();
854 Eigen::VectorXd distance(dimensions);
855
856 PetscInt colOwnerRangeCBegin, colOwnerRangeCEnd;
857 std::tie(colOwnerRangeCBegin, colOwnerRangeCEnd) = _matrixC.ownerRangeColumn();
858
859 VertexData vertexData(n - localPolyparams);
860
861 size_t local_row = 0;
862 // -- PREALLOCATES THE POLYNOMIAL PART OF THE MATRIX --
863 if (_polynomial == Polynomial::ON) {
864 for (local_row = 0; local_row < localPolyparams; local_row++) {
865 d_nnz[local_row] = colOwnerRangeCEnd - colOwnerRangeCBegin;
866 o_nnz[local_row] = _matrixC.getSize().first - d_nnz[local_row];
867 }
868 }
869
870 for (const mesh::Vertex &inVertex : inMesh->vertices()) {
871 if (not inVertex.isOwner())
872 continue;
873
874 PetscInt col = polyparams;
875 PetscInt const global_row = local_row + _matrixC.ownerRange().first;
876 d_nnz[local_row] = 0;
877 o_nnz[local_row] = 0;
878
879 // -- PREALLOCATES THE COEFFICIENTS --
880
881 for (auto const i : inMesh->index().getVerticesInsideBox(inVertex, supportRadius)) {
882 const mesh::Vertex &vj = inMesh->vertex(i);
883
884 PetscInt mappedCol = vj.getGlobalIndex() + polyparams;
885 AOApplicationToPetsc(_AOmapping, 1, &mappedCol); // likely not efficient in the inner loop
886
887 if (global_row > mappedCol) // Skip, since we are below the diagonal
888 continue;
889
890 distance = inVertex.getCoords() - vj.getCoords();
891 for (int d = 0; d < dimensions; d++)
892 if (this->_deadAxis[d])
893 distance[d] = 0;
894
895 double const norm = distance.norm();
896 if (supportRadius > norm or col == global_row) {
897 vertexData[local_row - localPolyparams].emplace_back(vj.getGlobalIndex() + polyparams, norm);
898 if (mappedCol >= colOwnerRangeCBegin and mappedCol < colOwnerRangeCEnd)
899 d_nnz[local_row]++;
900 else
901 o_nnz[local_row]++;
902 }
903 col++;
904 }
905 local_row++;
906 }
907
908 if (_commState->size() == 1) {
909 MatSeqSBAIJSetPreallocation(_matrixC, _matrixC.blockSize(), 0, d_nnz.data());
910 } else {
911 MatMPISBAIJSetPreallocation(_matrixC, _matrixC.blockSize(), 0, d_nnz.data(), 0, o_nnz.data());
912 }
913 MatSetOption(_matrixC, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);
914
915 ePreallocC.stop();
916
917 return vertexData;
918}
919
920template <typename RADIAL_BASIS_FUNCTION_T>
923{
924 PRECICE_INFO("Using tree-based preallocation for matrix A");
925 precice::profiling::Event ePreallocA("map.pet.preallocA.From" + this->input()->getName() + "To" + this->output()->getName(), profiling::Synchronize);
926
927 PetscInt ownerRangeABegin, ownerRangeAEnd, colOwnerRangeABegin, colOwnerRangeAEnd;
928 PetscInt const outputSize = _matrixA.getLocalSize().first;
929 double const supportRadius = this->_basisFunction.getSupportRadius();
930
931 std::tie(ownerRangeABegin, ownerRangeAEnd) = _matrixA.ownerRange();
932 std::tie(colOwnerRangeABegin, colOwnerRangeAEnd) = _matrixA.ownerRangeColumn();
933 const int dimensions = this->input()->getDimensions();
934
935 std::vector<PetscInt> d_nnz(outputSize), o_nnz(outputSize);
936 Eigen::VectorXd distance(dimensions);
937
938 // Contains localRow<localCols<colPosition, distance>>>
939 VertexData vertexData(outputSize);
940
941 for (int localRow = 0; localRow < ownerRangeAEnd - ownerRangeABegin; ++localRow) {
942 d_nnz[localRow] = 0;
943 o_nnz[localRow] = 0;
944 PetscInt col = 0;
945 mesh::Vertex const &oVertex = outMesh->vertex(localRow);
946
947 // -- PREALLOCATE THE POLYNOM PART OF THE MATRIX --
948 // col does not need mapping here, because the first polyparams col are always identity mapped
949 if (_polynomial == Polynomial::ON) {
950 if (col >= colOwnerRangeABegin and col < colOwnerRangeAEnd)
951 d_nnz[localRow]++;
952 else
953 o_nnz[localRow]++;
954 col++;
955
956 for (int dim = 0; dim < dimensions; dim++) {
957 if (not this->_deadAxis[dim]) {
958 if (col >= colOwnerRangeABegin and col < colOwnerRangeAEnd)
959 d_nnz[localRow]++;
960 else
961 o_nnz[localRow]++;
962 col++;
963 }
964 }
965 }
966
967 // -- PREALLOCATE THE COEFFICIENTS --
968 for (auto i : inMesh->index().getVerticesInsideBox(oVertex, supportRadius)) {
969 const mesh::Vertex &inVertex = inMesh->vertex(i);
970 distance = oVertex.getCoords() - inVertex.getCoords();
971
972 for (int d = 0; d < dimensions; d++)
973 if (this->_deadAxis[d])
974 distance[d] = 0;
975
976 double const norm = distance.norm();
977 if (supportRadius > norm) {
978 col = inVertex.getGlobalIndex() + polyparams;
979 vertexData[localRow].emplace_back(col, norm);
980
981 AOApplicationToPetsc(_AOmapping, 1, &col);
982 if (col >= colOwnerRangeABegin and col < colOwnerRangeAEnd)
983 d_nnz[localRow]++;
984 else
985 o_nnz[localRow]++;
986 }
987 }
988 }
989 if (_commState->size() == 1) {
990 MatSeqAIJSetPreallocation(_matrixA, 0, d_nnz.data());
991 } else {
992 MatMPIAIJSetPreallocation(_matrixA, 0, d_nnz.data(), 0, o_nnz.data());
993 }
994 MatSetOption(_matrixA, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);
995
996 ePreallocA.stop();
997 return vertexData;
998}
999
1000} // namespace mapping
1001} // namespace precice
1002
1003#endif // PRECICE_NO_PETSC
std::size_t distance
std::ostream & out
#define PRECICE_ERROR(...)
Definition LogMacros.hpp:15
#define PRECICE_WARN(...)
Definition LogMacros.hpp:11
#define PRECICE_DEBUG(...)
Definition LogMacros.hpp:64
#define PRECICE_DEBUG_IF(condition,...)
Definition LogMacros.hpp:66
#define PRECICE_TRACE(...)
Definition LogMacros.hpp:95
#define PRECICE_INFO(...)
Definition LogMacros.hpp:13
#define PRECICE_ASSERT(...)
Definition assertion.hpp:87
This class provides a lightweight logger.
Definition Logger.hpp:16
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 using the Petsc library to solve the resulting system.
size_t polyparams
Number of coefficients for the integrated polynomial. Depends on dimension and number of dead dimensi...
petsc::Matrix _matrixC
Interpolation system matrix. Evaluated basis function on the input mesh.
void mapConservative(const time::Sample &inData, Eigen::VectorXd &outData) final override
Maps data using a conservative constraint.
bool useRescaling
Toggles use of rescaled basis functions, only active when Polynomial == SEPARATE.
VertexData bgPreallocationMatrixC(const mesh::PtrMesh inMesh)
Preallocate matrix C and saves the coefficients using a boost::geometry spatial tree for neighbor sea...
void printMappingInfo(int dim) const
Prints an INFO about the current mapping.
void mapConsistent(const time::Sample &inData, Eigen::VectorXd &outData) final override
Maps data using a consistent constraint.
petsc::KSPSolver _solver
Used to solve matrixC for the RBF weighting factors.
PetRadialBasisFctMapping(Mapping::Constraint constraint, int dimensions, const RADIAL_BASIS_FUNCTION_T &function, std::array< bool, 3 > deadAxis, double solverRtol=1e-9, Polynomial polynomial=Polynomial::SEPARATE)
Constructor.
~PetRadialBasisFctMapping() override
Deletes the PETSc objects and the _deadAxis array.
void computeMapping() final override
Computes the mapping coefficients from the in- and output mesh.
std::string getName() const final override
name of the rbf mapping
void clear() final override
Removes a computed mapping.
VertexData bgPreallocationMatrixA(const mesh::PtrMesh inMesh, const mesh::PtrMesh outMesh)
petsc::Matrix _matrixQ
Vandermonde Matrix for linear polynomial, constructed from vertices of the input mesh.
void loadInitialGuessForDim(int dimension, int allDimensions, petsc::Vector &destination)
petsc::Matrix _matrixA
Interpolation evaluation matrix. Evaluated basis function on the output mesh.
petsc::KSPSolver _QRsolver
Used to solve the under-determined system for the separated polynomial.
void storeInitialGuessForDim(int dimension, int allDimensions, petsc::Vector &source)
Polynomial _polynomial
Toggles the use of the additional polynomial.
petsc::Matrix _matrixV
Coordinates of the output mesh to evaluate the separated polynomial.
utils::Parallel::CommStatePtr _commState
The CommState used to communicate.
size_t sepPolyparams
Number of coefficients for the separated polynomial. Depends on dimension and number of dead dimensio...
size_t localPolyparams
Equals to polyparams if rank == 0. Is 0 everywhere else.
petsc::Vector oneInterpolant
Interpolant of g(x) = 1 evaluated at output sites, used for rescaling.
int getPolynomialParameters() const
Computes the number of polynomial degrees of freedom based on the problem dimension and the dead axis...
Vertex of a mesh.
Definition Vertex.hpp:16
double coord(int index) const
Returns a coordinate of a vertex.
Definition Vertex.hpp:128
Eigen::VectorXd getCoords() const
Returns the coordinates of the vertex.
Definition Vertex.hpp:116
int getGlobalIndex() const
Globally unique index.
Definition Vertex.cpp:12
void stop()
Stops a running event.
Definition Event.cpp:37
void addData(std::string_view key, int value)
Adds named integer data, associated to an event.
Definition Event.cpp:48
Vector & copyTo(precice::span< double > destination)
Definition Petsc.cpp:379
PetscInt getLocalSize() const
Definition Petsc.cpp:313
Vector & copyFrom(precice::span< const double > source)
Definition Petsc.cpp:362
PetscInt getSize() const
Definition Petsc.cpp:304
T data(T... args)
T emplace_back(T... args)
T max(T... args)
Polynomial
How to handle the polynomial?
static constexpr SynchronizeTag Synchronize
Convenience instance of the SynchronizeTag.
Definition Event.hpp:21
PETSc related utilities.
Definition Petsc.cpp:121
void destroy(ISLocalToGlobalMapping *IS)
Destroys an ISLocalToGlobalMapping, if IS is not null and PetscIsInitialized.
Definition Petsc.cpp:796
Main namespace of the precice library.
T next(T... args)
T push_back(T... args)
T reserve(T... args)
T resize(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 tie(T... args)