preCICE v3.1.1
Loading...
Searching...
No Matches
MappingConfiguration.cpp
Go to the documentation of this file.
2#include <Eigen/Core>
3#include <algorithm>
4#include <cstring>
5#include <list>
6#include <memory>
7#include <ostream>
8#include <utility>
9#include <variant>
10#include "logging/LogMacros.hpp"
14#include "mapping/Mapping.hpp"
24#include "mesh/Mesh.hpp"
27#include "utils/Parallel.hpp"
28#include "utils/Petsc.hpp"
29#include "utils/assertion.hpp"
30#include "xml/ConfigParser.hpp"
31#include "xml/XMLAttribute.hpp"
32#include "xml/XMLTag.hpp"
33namespace precice::mapping {
34
35namespace {
36
37// given a list of subtags and parent tags, this function adds all subtags to all
38// parent tags
39void addSubtagsToParents(std::list<xml::XMLTag> &subtags,
41{
42 for (auto &p : parents) {
43 p.addSubtags(subtags);
44 }
45}
46
47// this function uses std::variant in order to add attributes of any type (double, string, bool)
48// to all tags in the list of tags \p storage.
49using variant_t = std::variant<xml::XMLAttribute<double>, xml::XMLAttribute<std::string>, xml::XMLAttribute<bool>, xml::XMLAttribute<int>>;
50template <typename TagStorage>
51void addAttributes(TagStorage &storage, const std::vector<variant_t> &attributes)
52{
53 for (auto &s : storage) {
54 for (auto &a : attributes)
55 std::visit([&s](auto &&arg) { s.addAttribute(arg); }, a);
56 }
57}
58
59// Enum required for the RBF instantiations
60enum struct RBFBackend {
61 Eigen,
62 PETSc,
63 Ginkgo,
64 PUM
65};
66
67// Helper in order to resolve the template instantiations.
68// Only the template specializations are of interest
69template <RBFBackend Backend, typename RBF>
70struct BackendSelector {
71 typedef RBF type;
72};
73
74// Specialization for the RBF Eigen backend
75template <typename RBF>
76struct BackendSelector<RBFBackend::Eigen, RBF> {
77 typedef mapping::RadialBasisFctMapping<RadialBasisFctSolver<RBF>> type;
78};
79
80// Specialization for the PETSc RBF backend
81#ifndef PRECICE_NO_PETSC
82template <typename RBF>
83struct BackendSelector<RBFBackend::PETSc, RBF> {
84 typedef mapping::PetRadialBasisFctMapping<RBF> type;
85};
86#endif
87
88// Specialization for the Ginkgo RBF backend
89#ifndef PRECICE_NO_GINKGO
90template <typename RBF>
91struct BackendSelector<RBFBackend::Ginkgo, RBF> {
92 typedef mapping::RadialBasisFctMapping<GinkgoRadialBasisFctSolver<RBF>, MappingConfiguration::GinkgoParameter> type;
93};
94#endif
95// Specialization for the RBF PUM backend
96template <typename RBF>
97struct BackendSelector<RBFBackend::PUM, RBF> {
98 typedef mapping::PartitionOfUnityMapping<RBF> type;
99};
100
101// Variant holding all available RBF classes
103
104// The actual instantiation of the mapping class, which is called by the visitor \ref getRBFMapping
105template <RBFBackend T, typename RADIAL_BASIS_FUNCTION_T, typename... Args>
106PtrMapping instantiateRBFMapping(mapping::Mapping::Constraint &constraint, int dimension, RADIAL_BASIS_FUNCTION_T function,
107 Args &&... args)
108{
109 return PtrMapping(new typename BackendSelector<T, RADIAL_BASIS_FUNCTION_T>::type(constraint, dimension, function, std::forward<Args>(args)...));
110}
111
112// Constructs the RBF function based on the functionType
113rbf_variant_t constructRBF(BasisFunction functionType, double supportRadius, double shapeParameter)
114{
115 switch (functionType) {
117 return mapping::CompactPolynomialC0(supportRadius);
118 }
120 return mapping::CompactPolynomialC2(supportRadius);
121 }
123 return mapping::CompactPolynomialC4(supportRadius);
124 }
126 return mapping::CompactPolynomialC6(supportRadius);
127 }
129 return mapping::CompactPolynomialC8(supportRadius);
130 }
132 return mapping::CompactThinPlateSplinesC2(supportRadius);
133 }
136 }
138 return mapping::VolumeSplines();
139 }
141 return mapping::Multiquadrics(shapeParameter);
142 }
144 return mapping::InverseMultiquadrics(shapeParameter);
145 }
147 return mapping::Gaussian(shapeParameter);
148 }
149 default:
150 PRECICE_UNREACHABLE("No instantiation was found for the selected basis function.");
151 }
152}
153
154// The actual instantion helper, which avoids enumerating all mapping implementations (more will come) with all RBF kernels
155// The first three arguments of the constructor are prescribed: constraint, dimension and the RBF function object, all other
156// constructor arguments are just forwareded. The first argument (BasisFunction) indicates then the actual instantiation to return.
157template <RBFBackend T, typename... Args>
158PtrMapping getRBFMapping(BasisFunction functionType, mapping::Mapping::Constraint &constraint, int dimension, double supportRadius, double shapeParameter,
159 Args &&... args)
160{
161 // First, construct the RBF function
162 auto functionVariant = constructRBF(functionType, supportRadius, shapeParameter);
163 // ... and instantiate the corresponding RBF mapping class
164 return std::visit([&](auto &&func) { return instantiateRBFMapping<T>(constraint, dimension, func, std::forward<Args>(args)...); }, functionVariant);
165}
166} // namespace
167
169 xml::XMLTag & parent,
170 mesh::PtrMeshConfiguration meshConfiguration)
171 : _meshConfig(std::move(meshConfiguration))
172{
174 using namespace xml;
175
176 // First, we create the available tags
177 XMLTag::Occurrence occ = XMLTag::OCCUR_ARBITRARY;
178 std::list<XMLTag> projectionTags{
179 XMLTag{*this, TYPE_NEAREST_NEIGHBOR, occ, TAG}.setDocumentation("Nearest-neighbour mapping which uses a rstar-spacial index tree to index meshes and run nearest-neighbour queries."),
180 XMLTag{*this, TYPE_NEAREST_PROJECTION, occ, TAG}.setDocumentation("Nearest-projection mapping which uses a rstar-spacial index tree to index meshes and locate the nearest projections."),
181 XMLTag{*this, TYPE_NEAREST_NEIGHBOR_GRADIENT, occ, TAG}.setDocumentation("Nearest-neighbor-gradient mapping which uses nearest-neighbor mapping with an additional linear approximation using gradient data."),
182 XMLTag{*this, TYPE_LINEAR_CELL_INTERPOLATION, occ, TAG}.setDocumentation("Linear cell interpolation mapping which uses a rstar-spacial index tree to index meshes and locate the nearest cell. Only supports 2D meshes.")};
183 std::list<XMLTag> rbfDirectTags{
184 XMLTag{*this, TYPE_RBF_GLOBAL_DIRECT, occ, TAG}.setDocumentation("Radial-basis-function mapping using a direct solver with a gather-scatter parallelism.")};
185 std::list<XMLTag> rbfIterativeTags{
186 XMLTag{*this, TYPE_RBF_GLOBAL_ITERATIVE, occ, TAG}.setDocumentation("Radial-basis-function mapping using an iterative solver with a distributed parallelism.")};
187 std::list<XMLTag> pumDirectTags{
188 XMLTag{*this, TYPE_RBF_PUM_DIRECT, occ, TAG}.setDocumentation("Radial-basis-function mapping using a partition of unity method, which supports a distributed parallelism.")};
189 std::list<XMLTag> rbfAliasTag{
190 XMLTag{*this, TYPE_RBF_ALIAS, occ, TAG}.setDocumentation("Alias tag, which auto-selects a radial-basis-function mapping depending on the simulation parameter,")};
191 std::list<XMLTag> geoMultiscaleTags{
192 XMLTag{*this, TYPE_AXIAL_GEOMETRIC_MULTISCALE, occ, TAG}.setDocumentation("Axial geometric multiscale mapping between one 1D and multiple 3D vertices."),
193 XMLTag{*this, TYPE_RADIAL_GEOMETRIC_MULTISCALE, occ, TAG}.setDocumentation("Radial geometric multiscale mapping between multiple 1D and multiple 3D vertices, distributed along a principle axis.")};
194
195 // List of all attributes with corresponding documentation
196 auto attrDirection = XMLAttribute<std::string>(ATTR_DIRECTION)
197 .setOptions({DIRECTION_WRITE, DIRECTION_READ})
198 .setDocumentation("Write mappings map written data prior to communication, thus in the same participant who writes the data. "
199 "Read mappings map received data after communication, thus in the same participant who reads the data.");
200
201 auto attrFromMesh = XMLAttribute<std::string>(ATTR_FROM)
202 .setDocumentation("The mesh to map the data from.");
203
204 auto attrToMesh = XMLAttribute<std::string>(ATTR_TO)
205 .setDocumentation("The mesh to map the data to.");
206
207 auto attrConstraint = XMLAttribute<std::string>(ATTR_CONSTRAINT)
208 .setDocumentation("Use conservative to conserve the nodal sum of the data over the interface (needed e.g. for force mapping). Use consistent for normalized quantities such as temperature or pressure. Use scaled-consistent-surface or scaled-consistent-volume for normalized quantities where conservation of integral values (surface or volume) is needed (e.g. velocities when the mass flow rate needs to be conserved). Mesh connectivity is required to use scaled-consistent.")
210 auto attrXDead = makeXMLAttribute(ATTR_X_DEAD, false)
211 .setDocumentation("If set to true, the x axis will be ignored for the mapping");
212 auto attrYDead = makeXMLAttribute(ATTR_Y_DEAD, false)
213 .setDocumentation("If set to true, the y axis will be ignored for the mapping");
214 auto attrZDead = makeXMLAttribute(ATTR_Z_DEAD, false)
215 .setDocumentation("If set to true, the z axis will be ignored for the mapping");
216 auto attrPolynomial = makeXMLAttribute(ATTR_POLYNOMIAL, POLYNOMIAL_SEPARATE)
217 .setDocumentation("Toggles use of the global polynomial")
219 auto attrPumPolynomial = makeXMLAttribute(ATTR_POLYNOMIAL, POLYNOMIAL_SEPARATE)
220 .setDocumentation("Toggles use a local (per cluster) polynomial")
222
223 auto attrSolverRtol = makeXMLAttribute(ATTR_SOLVER_RTOL, 1e-9)
224 .setDocumentation("Solver relative tolerance for convergence");
225 // TODO: Discuss whether we wanto to introduce this attribute
226 // auto attrMaxIterations = makeXMLAttribute(ATTR_MAX_ITERATIONS, 1e6)
227 // .setDocumentation("Maximum number of iterations of the solver");
228
229 auto verticesPerCluster = XMLAttribute<int>(ATTR_VERTICES_PER_CLUSTER, 50)
230 .setDocumentation("Average number of vertices per cluster (partition) applied in the rbf partition of unity method.");
231 auto relativeOverlap = makeXMLAttribute(ATTR_RELATIVE_OVERLAP, 0.15)
232 .setDocumentation("Value between 0 and 1 indicating the relative overlap between clusters. A value of 0.15 is usually a good trade-off between accuracy and efficiency.");
233 auto projectToInput = XMLAttribute<bool>(ATTR_PROJECT_TO_INPUT, true)
234 .setDocumentation("If enabled, places the cluster centers at the closest vertex of the input mesh. Should be enabled in case of non-uniform point distributions such as for shell structures.");
235
236 auto attrGeoMultiscaleType = XMLAttribute<std::string>(ATTR_GEOMETRIC_MULTISCALE_TYPE)
237 .setDocumentation("Type of geometric multiscale mapping. Either 'spread' or 'collect'.")
239 auto attrGeoMultiscaleAxis = XMLAttribute<std::string>(ATTR_GEOMETRIC_MULTISCALE_AXIS)
240 .setDocumentation("Principle axis along which geometric multiscale mapping is performed.")
242 auto attrGeoMultiscaleRadius = XMLAttribute<double>(ATTR_GEOMETRIC_MULTISCALE_RADIUS)
243 .setDocumentation("Radius of the circular interface between the 1D and 3D participant.");
244
245 // Add the relevant attributes to the relevant tags
246 addAttributes(projectionTags, {attrFromMesh, attrToMesh, attrDirection, attrConstraint});
247 addAttributes(rbfDirectTags, {attrFromMesh, attrToMesh, attrDirection, attrConstraint, attrPolynomial, attrXDead, attrYDead, attrZDead});
248 addAttributes(rbfIterativeTags, {attrFromMesh, attrToMesh, attrDirection, attrConstraint, attrPolynomial, attrXDead, attrYDead, attrZDead, attrSolverRtol});
249 addAttributes(pumDirectTags, {attrFromMesh, attrToMesh, attrDirection, attrConstraint, attrPumPolynomial, verticesPerCluster, relativeOverlap, projectToInput});
250 addAttributes(rbfAliasTag, {attrFromMesh, attrToMesh, attrDirection, attrConstraint, attrXDead, attrYDead, attrZDead});
251 addAttributes(geoMultiscaleTags, {attrFromMesh, attrToMesh, attrDirection, attrConstraint, attrGeoMultiscaleType, attrGeoMultiscaleAxis, attrGeoMultiscaleRadius});
252
253 // Now we take care of the subtag executor. We repeat some of the subtags in order to add individual documentation
254 XMLTag::Occurrence once = XMLTag::OCCUR_NOT_OR_ONCE;
255 // TODO, make type an int
256 auto attrDeviceId = makeXMLAttribute(ATTR_DEVICE_ID, static_cast<int>(0))
257 .setDocumentation("Specifies the ID of the GPU that should be used for the Ginkgo GPU backend.");
258 auto attrNThreads = makeXMLAttribute(ATTR_N_THREADS, static_cast<int>(0))
259 .setDocumentation("Specifies the number of threads for the OpenMP executor that should be used for the Ginkgo OpenMP backend. If a value of \"0\" is set, preCICE doesn't set the number of threads and the default behavior of OpenMP applies.");
260
261 // First, we have the executors for the direct solvers
262 {
263 std::list<XMLTag> cpuExecutor{
264 XMLTag{*this, EXECUTOR_CPU, once, SUBTAG_EXECUTOR}.setDocumentation("The default executor, which uses a single-core CPU with a gather-scatter parallelism.")};
265 std::list<XMLTag> deviceExecutors{
266 XMLTag{*this, EXECUTOR_CUDA, once, SUBTAG_EXECUTOR}.setDocumentation("Cuda (Nvidia) executor, which uses cuSolver/Ginkgo and a direct QR decomposition with a gather-scatter parallelism."),
267 XMLTag{*this, EXECUTOR_HIP, once, SUBTAG_EXECUTOR}.setDocumentation("Hip (AMD/Nvidia) executor, which uses hipSolver/Ginkgo and a direct QR decomposition with a gather-scatter parallelism.")};
268
269 addAttributes(deviceExecutors, {attrDeviceId});
270 addSubtagsToParents(cpuExecutor, rbfDirectTags);
271 addSubtagsToParents(deviceExecutors, rbfDirectTags);
272 }
273 // Second, the executors for the iterative solver
274 {
275 std::list<XMLTag> cpuExecutor{
276 XMLTag{*this, EXECUTOR_CPU, once, SUBTAG_EXECUTOR}.setDocumentation("The default executor relying on PETSc, which uses CPUs and distributed memory parallelism via MPI.")};
277 std::list<XMLTag> deviceExecutors{
278 XMLTag{*this, EXECUTOR_CUDA, once, SUBTAG_EXECUTOR}.setDocumentation("Cuda (Nvidia) executor, which uses Ginkgo with a gather-scatter parallelism."),
279 XMLTag{*this, EXECUTOR_HIP, once, SUBTAG_EXECUTOR}.setDocumentation("Hip (AMD/Nvidia) executor, which uses hipSolver with a gather-scatter parallelism.")};
280 std::list<XMLTag> ompExecutor{
281 XMLTag{*this, EXECUTOR_OMP, once, SUBTAG_EXECUTOR}.setDocumentation("OpenMP executor, which uses Ginkgo with a gather-scatter parallelism.")};
282
283 addAttributes(deviceExecutors, {attrDeviceId});
284 addAttributes(ompExecutor, {attrNThreads});
285 addSubtagsToParents(cpuExecutor, rbfIterativeTags);
286 addSubtagsToParents(deviceExecutors, rbfIterativeTags);
287 addSubtagsToParents(ompExecutor, rbfIterativeTags);
288 }
289 {
290 std::list<XMLTag> cpuExecutor{
291 XMLTag{*this, EXECUTOR_CPU, once, SUBTAG_EXECUTOR}.setDocumentation("The default (and currently only) executor using a CPU and a distributed memory parallelism via MPI.")};
292 addSubtagsToParents(cpuExecutor, pumDirectTags);
293 }
294 // The alias tag doesn't receive the subtag at all
295
296 // Now we take care of the subtag basis function
297 // First, we have the tags using a support radius
298 std::list<XMLTag> supportRadiusRBF{
299 XMLTag{*this, RBF_CPOLYNOMIAL_C0, once, SUBTAG_BASIS_FUNCTION}.setDocumentation("Wendland C0 function"),
300 XMLTag{*this, RBF_CPOLYNOMIAL_C2, once, SUBTAG_BASIS_FUNCTION}.setDocumentation("Wendland C2 function"),
301 XMLTag{*this, RBF_CPOLYNOMIAL_C4, once, SUBTAG_BASIS_FUNCTION}.setDocumentation("Wendland C4 function"),
302 XMLTag{*this, RBF_CPOLYNOMIAL_C6, once, SUBTAG_BASIS_FUNCTION}.setDocumentation("Wendland C6 function"),
303 XMLTag{*this, RBF_CPOLYNOMIAL_C8, once, SUBTAG_BASIS_FUNCTION}.setDocumentation("Wendland C8 function"),
304 XMLTag{*this, RBF_CTPS_C2, once, SUBTAG_BASIS_FUNCTION}.setDocumentation("Compact thin-plate-spline C2")};
305
306 auto attrSupportRadius = XMLAttribute<double>(ATTR_SUPPORT_RADIUS)
307 .setDocumentation("Support radius of each RBF basis function (global choice).");
308
309 addAttributes(supportRadiusRBF, {attrSupportRadius});
310 addSubtagsToParents(supportRadiusRBF, rbfIterativeTags);
311 addSubtagsToParents(supportRadiusRBF, rbfDirectTags);
312 addSubtagsToParents(supportRadiusRBF, pumDirectTags);
313 addSubtagsToParents(supportRadiusRBF, rbfAliasTag);
314
315 // Now the tags using a shape parameter
316 std::list<XMLTag> shapeParameterRBF{
317 XMLTag{*this, RBF_MULTIQUADRICS, once, SUBTAG_BASIS_FUNCTION}.setDocumentation("Multiquadrics"),
318 XMLTag{*this, RBF_INV_MULTIQUADRICS, once, SUBTAG_BASIS_FUNCTION}.setDocumentation("Inverse multiquadrics")};
319
320 auto attrShapeParam = XMLAttribute<double>(ATTR_SHAPE_PARAM)
321 .setDocumentation("Specific shape parameter for RBF basis function.");
322
323 addAttributes(shapeParameterRBF, {attrShapeParam});
324 addSubtagsToParents(shapeParameterRBF, rbfIterativeTags);
325 addSubtagsToParents(shapeParameterRBF, rbfDirectTags);
326 addSubtagsToParents(shapeParameterRBF, pumDirectTags);
327 addSubtagsToParents(shapeParameterRBF, rbfAliasTag);
328
329 // For the Gaussian, we need default values as the user can pass a support radius or a shape parameter
330 std::list<XMLTag> GaussRBF{
331 XMLTag{*this, RBF_GAUSSIAN, once, SUBTAG_BASIS_FUNCTION}.setDocumentation("Gaussian basis function accepting a support radius or a shape parameter.")};
332 attrShapeParam.setDefaultValue(std::numeric_limits<double>::quiet_NaN());
333 attrSupportRadius.setDefaultValue(std::numeric_limits<double>::quiet_NaN());
334 addAttributes(GaussRBF, {attrShapeParam, attrSupportRadius});
335 addSubtagsToParents(GaussRBF, rbfIterativeTags);
336 addSubtagsToParents(GaussRBF, rbfDirectTags);
337 addSubtagsToParents(GaussRBF, pumDirectTags);
338 addSubtagsToParents(GaussRBF, rbfAliasTag);
339
340 // tags without an attribute
341 std::list<XMLTag> attributelessRBFs{
342 XMLTag{*this, RBF_TPS, once, SUBTAG_BASIS_FUNCTION}.setDocumentation("Thin-plate-splines"),
343 XMLTag{*this, RBF_VOLUME_SPLINES, once, SUBTAG_BASIS_FUNCTION}.setDocumentation("Volume splines")};
344
345 addSubtagsToParents(attributelessRBFs, rbfIterativeTags);
346 addSubtagsToParents(attributelessRBFs, rbfDirectTags);
347 addSubtagsToParents(attributelessRBFs, pumDirectTags);
348 addSubtagsToParents(attributelessRBFs, rbfAliasTag);
349
350 // Add all tags to the mapping tag
351 parent.addSubtags(projectionTags);
352 parent.addSubtags(rbfIterativeTags);
353 parent.addSubtags(rbfDirectTags);
354 parent.addSubtags(pumDirectTags);
355 parent.addSubtags(rbfAliasTag);
356 parent.addSubtags(geoMultiscaleTags);
357}
358
360 bool experimental)
361{
362 _experimental = experimental;
363}
364
366 const xml::ConfigurationContext &context,
367 xml::XMLTag & tag)
368{
369 PRECICE_TRACE(tag.getName());
370 if (tag.getNamespace() == TAG) {
371 // Mandatory tags
375 std::string type = tag.getName();
377
378 // optional tags
379 // We set here default values, but their actual value doesn't really matter.
380 // It's just for the mapping methods, which do not use these attributes at all.
381 bool xDead = tag.getBooleanAttributeValue(ATTR_X_DEAD, false);
382 bool yDead = tag.getBooleanAttributeValue(ATTR_Y_DEAD, false);
383 bool zDead = tag.getBooleanAttributeValue(ATTR_Z_DEAD, false);
384 double solverRtol = tag.getDoubleAttributeValue(ATTR_SOLVER_RTOL, 1e-9);
386
387 // geometric multiscale related tags
390 double multiscaleRadius = tag.getDoubleAttributeValue(ATTR_GEOMETRIC_MULTISCALE_RADIUS, 1.0);
391
393 PRECICE_CHECK(_experimental, "Axial geometric multiscale is experimental and the configuration can change between minor releases. Set experimental=\"on\" in the precice-configuration tag.");
394 }
395
396 if (type == TYPE_AXIAL_GEOMETRIC_MULTISCALE && context.size > 1) {
397 throw std::runtime_error{"Axial geometric multiscale mapping is not available for parallel participants."};
398 }
399
400 if (type == TYPE_RADIAL_GEOMETRIC_MULTISCALE && context.size > 1) {
401 throw std::runtime_error{"Radial geometric multiscale mapping is not available for parallel participants."};
402 }
403
404 // pum related tags
405 int verticesPerCluster = tag.getIntAttributeValue(ATTR_VERTICES_PER_CLUSTER, 100);
406 double relativeOverlap = tag.getDoubleAttributeValue(ATTR_RELATIVE_OVERLAP, 0.3);
407 bool projectToInput = tag.getBooleanAttributeValue(ATTR_PROJECT_TO_INPUT, true);
408
409 // Convert raw string into enum types as the constructors take enums
410 if (constraint == CONSTRAINT_CONSERVATIVE) {
412 } else if (constraint == CONSTRAINT_CONSISTENT) {
414 } else if (constraint == CONSTRAINT_SCALED_CONSISTENT_SURFACE) {
416 } else if (constraint == CONSTRAINT_SCALED_CONSISTENT_VOLUME) {
418 } else {
419 PRECICE_UNREACHABLE("Unknown mapping constraint \"{}\".", constraint);
420 }
421
422 ConfiguredMapping configuredMapping = createMapping(dir, type, fromMesh, toMesh, geoMultiscaleType, geoMultiscaleAxis, multiscaleRadius);
423
424 _rbfConfig = configureRBFMapping(type, strPolynomial, xDead, yDead, zDead, solverRtol, verticesPerCluster, relativeOverlap, projectToInput);
425
426 checkDuplicates(configuredMapping);
427 _mappings.push_back(configuredMapping);
428 } else if (tag.getNamespace() == SUBTAG_BASIS_FUNCTION) {
429
430 PRECICE_ASSERT(!_mappings.empty());
431 PRECICE_CHECK(_mappings.back().requiresBasisFunction == true, "A basis-function was defined for the mapping "
432 "from mesh \"{}\" to mesh \"{}\", but no basis-function is required for this mapping type. "
433 "Please remove the basis-function tag or configure an rbf mapping, which requires a basis-function.",
434 _mappings.back().fromMesh->getName(), _mappings.back().toMesh->getName());
435 // We can only set one subtag
436 PRECICE_CHECK(_rbfConfig.basisFunctionDefined == false, "More than one basis-function was defined for the mapping "
437 "from mesh \"{}\" to mesh \"{}\".",
438 _mappings.back().fromMesh->getName(), _mappings.back().toMesh->getName());
439
440 std::string basisFctName = tag.getName();
441 double supportRadius = tag.getDoubleAttributeValue(ATTR_SUPPORT_RADIUS, 0.);
442 double shapeParameter = tag.getDoubleAttributeValue(ATTR_SHAPE_PARAM, 0.);
443
446 // The Gaussian RBF is always treated as a shape-parameter RBF. Hence, we have to convert the support radius, if necessary
448 const bool exactlyOneSet = (std::isfinite(supportRadius) && !std::isfinite(shapeParameter)) ||
449 (std::isfinite(shapeParameter) && !std::isfinite(supportRadius));
450 PRECICE_CHECK(exactlyOneSet, "The specified parameters for the Gaussian RBF mapping are invalid. Please specify either a \"shape-parameter\" or a \"support-radius\".");
451
452 if (std::isfinite(supportRadius) && !std::isfinite(shapeParameter)) {
453 shapeParameter = std::sqrt(-std::log(Gaussian::cutoffThreshold)) / supportRadius;
454 }
455 }
456
457 _rbfConfig.supportRadius = supportRadius;
458 _rbfConfig.shapeParameter = shapeParameter;
459 } else if (tag.getNamespace() == SUBTAG_EXECUTOR) {
460 _executorConfig = std::make_unique<ExecutorConfiguration>();
461
462 if (tag.getName() == EXECUTOR_CPU) {
464 } else if (tag.getName() == EXECUTOR_CUDA) {
466 } else if (tag.getName() == EXECUTOR_HIP) {
468 } else if (tag.getName() == EXECUTOR_OMP) {
470 }
471
474 }
475}
476
478 const std::string &polynomial,
479 bool xDead, bool yDead, bool zDead,
480 double solverRtol,
481 double verticesPerCluster,
482 double relativeOverlap,
483 bool projectToInput) const
484{
486
487 if (type == TYPE_RBF_GLOBAL_ITERATIVE)
489 else if (type == TYPE_RBF_GLOBAL_DIRECT)
491 else if (type == TYPE_RBF_PUM_DIRECT)
493 else {
494 // Rather simple auto-selection (for now), consisting of the PUM Eigen backend
496
497 // A more sophisticated criterion here could take the globalNumberOfVertices into account
498 // (the mesh pointer is stored in the configuredMapping anyway), but this quantity is not yet computed
499 // during the configuration time.
500 }
501
502 if (polynomial == POLYNOMIAL_SEPARATE)
504 else if (polynomial == POLYNOMIAL_ON)
506 else if (polynomial == POLYNOMIAL_OFF)
508 else
509 PRECICE_UNREACHABLE("Unknown polynomial configuration.");
510
511 rbfConfig.deadAxis = {{xDead, yDead, zDead}};
512 rbfConfig.solverRtol = solverRtol;
513
514 rbfConfig.verticesPerCluster = verticesPerCluster;
515 rbfConfig.relativeOverlap = relativeOverlap;
516 rbfConfig.projectToInput = projectToInput;
517
518 return rbfConfig;
519}
520
522 const std::string &direction,
523 const std::string &type,
524 const std::string &fromMeshName,
525 const std::string &toMeshName,
526 const std::string &geoMultiscaleType,
527 const std::string &geoMultiscaleAxis,
528 const double & multiscaleRadius) const
529{
530 PRECICE_TRACE(direction, type);
531
532 ConfiguredMapping configuredMapping;
533 mesh::PtrMesh fromMesh(_meshConfig->getMesh(fromMeshName));
534 mesh::PtrMesh toMesh(_meshConfig->getMesh(toMeshName));
535 PRECICE_CHECK(fromMesh.get() != nullptr,
536 "Mesh \"{0}\" was not found while creating a mapping. "
537 "Please correct the from=\"{0}\" attribute.",
538 fromMeshName);
539 PRECICE_CHECK(toMesh.get() != nullptr,
540 "Mesh \"{0}\" was not found while creating a mapping. "
541 "Please correct the to=\"{0}\" attribute.",
542 toMeshName);
543
544 // Check for compatible mesh dimensions
545 PRECICE_CHECK(fromMesh->getDimensions() == toMesh->getDimensions(),
546 "Mapping between meshes of different dimensions is not allowed yet. "
547 "Please set the same dimensions attribute to meshes \"{}\" and \"{}\", "
548 "or choose different meshes.",
549 fromMeshName, toMeshName);
550
551 configuredMapping.fromMesh = fromMesh;
552 configuredMapping.toMesh = toMesh;
553
554 if (direction == DIRECTION_WRITE) {
555 configuredMapping.direction = WRITE;
556 } else if (direction == DIRECTION_READ) {
557 configuredMapping.direction = READ;
558 } else {
559 PRECICE_UNREACHABLE("Unknown mapping direction type \"{}\".", direction);
560 }
561
562 // Create all projection based mappings
563 if (type == TYPE_NEAREST_NEIGHBOR) {
564 configuredMapping.mapping = PtrMapping(new NearestNeighborMapping(constraintValue, fromMesh->getDimensions()));
565 } else if (type == TYPE_NEAREST_PROJECTION) {
566 configuredMapping.mapping = PtrMapping(new NearestProjectionMapping(constraintValue, fromMesh->getDimensions()));
567 } else if (type == TYPE_LINEAR_CELL_INTERPOLATION) {
568 configuredMapping.mapping = PtrMapping(new LinearCellInterpolationMapping(constraintValue, fromMesh->getDimensions()));
569 } else if (type == TYPE_NEAREST_NEIGHBOR_GRADIENT) {
570
571 // NNG is not applicable with the conservative constraint
573 "Nearest-neighbor-gradient mapping is not implemented using a \"conservative\" constraint. "
574 "Please select constraint=\" consistent\" or a different mapping method.");
575
576 configuredMapping.mapping = PtrMapping(new NearestNeighborGradientMapping(constraintValue, fromMesh->getDimensions()));
577
578 } else if (type == TYPE_AXIAL_GEOMETRIC_MULTISCALE) {
579
580 // Axial geometric multiscale is not applicable with the conservative constraint
582 "Axial geometric multiscale mapping is not implemented for the \"conservative\" constraint. "
583 "Please select constraint=\" consistent\" or a different mapping method.");
584
585 // Convert strings into enums
587 if (geoMultiscaleAxis == "x") {
589 } else if (geoMultiscaleAxis == "y") {
591 } else if (geoMultiscaleAxis == "z") {
593 } else {
594 PRECICE_UNREACHABLE("Unknown geometric multiscale axis \"{}\".", geoMultiscaleAxis);
595 }
596
598 if (geoMultiscaleType == "spread") {
600 } else if (geoMultiscaleType == "collect") {
602 } else {
603 PRECICE_UNREACHABLE("Unknown geometric multiscale type \"{}\".", geoMultiscaleType);
604 }
605
606 configuredMapping.mapping = PtrMapping(new AxialGeoMultiscaleMapping(constraintValue, fromMesh->getDimensions(), multiscaleType, multiscaleAxis, multiscaleRadius));
607
608 } else if (type == TYPE_RADIAL_GEOMETRIC_MULTISCALE) {
609
610 // Radial geometric multiscale is not applicable with the conservative constraint
612 "Radial geometric multiscale mapping is not implemented for the \"conservative\" constraint. "
613 "Please select constraint=\" consistent\" or a different mapping method.");
614
615 // Convert strings into enums
617 if (geoMultiscaleAxis == "x") {
619 } else if (geoMultiscaleAxis == "y") {
621 } else if (geoMultiscaleAxis == "z") {
623 } else {
624 PRECICE_UNREACHABLE("Unknown geometric multiscale axis \"{}\".", geoMultiscaleAxis);
625 }
626
628 if (geoMultiscaleType == "spread") {
630 } else if (geoMultiscaleType == "collect") {
632 } else {
633 PRECICE_UNREACHABLE("Unknown geometric multiscale type \"{}\".", geoMultiscaleType);
634 }
635
636 configuredMapping.mapping = PtrMapping(new RadialGeoMultiscaleMapping(constraintValue, fromMesh->getDimensions(), multiscaleType, multiscaleAxis));
637
638 } else {
639 // We need knowledge about the basis function in order to instantiate the rbf related mapping
641 configuredMapping.mapping = nullptr;
642 configuredMapping.configuredWithAliasTag = type == TYPE_RBF_ALIAS;
643 }
644
645 configuredMapping.requiresBasisFunction = requiresBasisFunction(type);
646
647 return configuredMapping;
648}
649
651{
652 for (const ConfiguredMapping &configuredMapping : _mappings) {
653 bool sameToMesh = mapping.toMesh->getName() == configuredMapping.toMesh->getName();
654 bool sameFromMesh = mapping.fromMesh->getName() == configuredMapping.fromMesh->getName();
655 bool sameMapping = sameToMesh && sameFromMesh;
656 PRECICE_CHECK(!sameMapping,
657 "There cannot be two mappings from mesh \"{}\" to mesh \"{}\". "
658 "Please remove one of the duplicated meshes. ",
659 mapping.fromMesh->getName(), mapping.toMesh->getName());
660 }
661}
662
664{
665 if (tag.getNamespace() == TAG) {
666 if (requiresBasisFunction(tag.getName())) {
667 PRECICE_CHECK(_rbfConfig.basisFunctionDefined, "No basis-function was defined for the \"{}\" mapping from mesh \"{}\" to mesh \"{}\".", tag.getName(), _mappings.back().fromMesh->getName(), _mappings.back().toMesh->getName());
668 if (!_executorConfig) {
669 _executorConfig = std::make_unique<ExecutorConfiguration>();
670 }
672 _executorConfig.reset();
673 }
674 PRECICE_ASSERT(_mappings.back().mapping != nullptr);
675 }
676}
677
679{
681 ConfiguredMapping &mapping = _mappings.back();
682 // Instantiate the RBF mapping classes
683 // We first categorize according to the executor
684 // 1. the CPU executor
687 mapping.mapping = getRBFMapping<RBFBackend::Eigen>(_rbfConfig.basisFunction, constraintValue, mapping.fromMesh->getDimensions(), _rbfConfig.supportRadius, _rbfConfig.shapeParameter, _rbfConfig.deadAxis, _rbfConfig.polynomial);
689#ifndef PRECICE_NO_PETSC
690 // for petsc initialization
692
694#else
695 PRECICE_CHECK(false, "The global-iterative RBF solver on a CPU requires a preCICE build with PETSc enabled.");
696#endif
699 } else {
700 PRECICE_UNREACHABLE("Unknown RBF solver.");
701 }
702 // 2. any other executor is configured via Ginkgo
703 } else {
704#ifndef PRECICE_NO_GINKGO
709 _ginkgoParameter.executor = "cuda-executor";
710#ifndef PRECICE_WITH_CUDA
711 PRECICE_CHECK(false, "The cuda-executor (configured for the mapping from mesh {} to mesh {}) requires a Ginkgo and preCICE build with Cuda enabled.", mapping.fromMesh->getName(), mapping.toMesh->getName());
712#endif
714 _ginkgoParameter.executor = "hip-executor";
715#ifndef PRECICE_WITH_HIP
716 PRECICE_CHECK(false, "The hip-executor (configured for the mapping from mesh {} to mesh {}) requires a Ginkgo and preCICE build with HIP enabled.", mapping.fromMesh->getName(), mapping.toMesh->getName());
717#endif
719 _ginkgoParameter.executor = "omp-executor";
721#ifndef PRECICE_WITH_OMP
722 PRECICE_CHECK(false, "The omp-executor (configured for the mapping from mesh {} to mesh {}) requires a Ginkgo and preCICE build with OpenMP enabled.", mapping.fromMesh->getName(), mapping.toMesh->getName());
723#endif
724 }
726 _ginkgoParameter.solver = "qr-solver";
728 _ginkgoParameter.solver = "cg-solver";
730 } else {
731 PRECICE_UNREACHABLE("Unknown solver type.");
732 }
734#else
735 PRECICE_CHECK(false, "The selected executor for the mapping from mesh {} to mesh {} requires a preCICE build with Ginkgo enabled.", mapping.fromMesh->getName(), mapping.toMesh->getName());
736#endif
737 }
738}
739
744
746{
747 return mappingType == TYPE_RBF_PUM_DIRECT || mappingType == TYPE_RBF_GLOBAL_DIRECT || mappingType == TYPE_RBF_GLOBAL_ITERATIVE || mappingType == TYPE_RBF_ALIAS;
748}
749
751{
752 BasisFunction basisFunction{};
753 if (basisFctName == RBF_TPS)
754 basisFunction = BasisFunction::ThinPlateSplines;
755 else if (basisFctName == RBF_MULTIQUADRICS)
756 basisFunction = BasisFunction::Multiquadrics;
757 else if (basisFctName == RBF_INV_MULTIQUADRICS)
759 else if (basisFctName == RBF_VOLUME_SPLINES)
760 basisFunction = BasisFunction::VolumeSplines;
761 else if (basisFctName == RBF_GAUSSIAN)
762 basisFunction = BasisFunction::Gaussian;
763 else if (basisFctName == RBF_CTPS_C2)
765 else if (basisFctName == RBF_CPOLYNOMIAL_C0)
766 basisFunction = BasisFunction::WendlandC0;
767 else if (basisFctName == RBF_CPOLYNOMIAL_C2)
768 basisFunction = BasisFunction::WendlandC2;
769 else if (basisFctName == RBF_CPOLYNOMIAL_C4)
770 basisFunction = BasisFunction::WendlandC4;
771 else if (basisFctName == RBF_CPOLYNOMIAL_C6)
772 basisFunction = BasisFunction::WendlandC6;
773 else if (basisFctName == RBF_CPOLYNOMIAL_C8)
774 basisFunction = BasisFunction::WendlandC8;
775 else
776 PRECICE_UNREACHABLE("Unknown basis function \"{}\".", basisFctName);
777 return basisFunction;
778}
779} // namespace precice::mapping
#define PRECICE_TRACE(...)
Definition LogMacros.hpp:95
#define PRECICE_CHECK(check,...)
Definition LogMacros.hpp:35
#define PRECICE_ASSERT(...)
Definition assertion.hpp:87
#define PRECICE_UNREACHABLE(...)
Definition assertion.hpp:95
Geometric multiscale mapping in axial direction.
MultiscaleType
Geometric multiscale nature of the mapping (spread or collect).
static constexpr double cutoffThreshold
Below that value the function is supposed to be zero. Defines the support radius if not explicitly gi...
Mapping using orthogonal projection to nearest triangle/edge/vertex and linear interpolation from pro...
virtual void xmlEndTagCallback(const xml::ConfigurationContext &context, xml::XMLTag &callingTag)
Callback function required for use of automatic configuration.
RBFConfiguration configureRBFMapping(const std::string &type, const std::string &polynomial, bool xDead, bool yDead, bool zDead, double solverRtol, double verticesPerCluster, double relativeOverlap, bool projectToInput) const
virtual void xmlTagCallback(const xml::ConfigurationContext &context, xml::XMLTag &callingTag)
Callback function required for use of automatic configuration.
const std::vector< ConfiguredMapping > & mappings()
Returns all configured mappings.
std::unique_ptr< ExecutorConfiguration > _executorConfig
bool requiresBasisFunction(const std::string &mappingType) const
BasisFunction parseBasisFunctions(const std::string &basisFctName) const
Given a basis function name (as a string), transforms the string into an enum of the BasisFunction.
std::vector< ConfiguredMapping > _mappings
ConfiguredMapping createMapping(const std::string &direction, const std::string &type, const std::string &fromMeshName, const std::string &toMeshName, const std::string &geoMultiscaleType, const std::string &geoMultiscaleAxis, const double &multiscaleRadius) const
void checkDuplicates(const ConfiguredMapping &mapping)
Check whether a mapping to and from the same mesh already exists.
const RBFConfiguration & rbfConfig() const
MappingConfiguration(xml::XMLTag &parent, mesh::PtrMeshConfiguration meshConfiguration)
Constraint
Specifies additional constraints for a mapping.
Definition Mapping.hpp:29
Mapping using nearest neighboring vertices and their local gradient values.
Mapping using nearest neighboring vertices.
Mapping using orthogonal projection to nearest triangle/edge/vertex and linear interpolation from pro...
Geometric multiscale mapping in radial direction.
MultiscaleType
Geometric multiscale type of the mapping (spread or collect).
static CommStatePtr current()
Returns an owning pointer to the current CommState.
Definition Parallel.cpp:147
static void initialize(utils::Parallel::Communicator comm)
Initializes the Petsc environment.
Definition Petsc.cpp:78
XMLAttribute & setOptions(std::vector< ATTRIBUTE_T > options)
const std::string & getName() const
XMLAttribute & setDocumentation(std::string documentation)
Sets a documentation string for the attribute.
Represents an XML tag to be configured automatically.
Definition XMLTag.hpp:31
const std::string & getNamespace() const
Returns xml namespace.
Definition XMLTag.hpp:159
void addSubtags(const Container &subtags)
Definition XMLTag.hpp:143
std::string getStringAttributeValue(const std::string &name, std::optional< std::string > default_value=std::nullopt) const
Definition XMLTag.cpp:142
bool getBooleanAttributeValue(const std::string &name, std::optional< bool > default_value=std::nullopt) const
Definition XMLTag.cpp:155
const std::string & getName() const
Returns name (without namespace).
Definition XMLTag.hpp:153
int getIntAttributeValue(const std::string &name, std::optional< int > default_value=std::nullopt) const
Definition XMLTag.cpp:129
double getDoubleAttributeValue(const std::string &name, std::optional< double > default_value=std::nullopt) const
Definition XMLTag.cpp:116
T get(T... args)
T isfinite(T... args)
T log(T... args)
contains data mapping from points to meshes.
std::shared_ptr< Mapping > PtrMapping
STL namespace.
T sqrt(T... args)
Direction direction
Direction of mapping (important to set input and output mesh).
bool configuredWithAliasTag
used the automatic rbf alias tag in order to set the mapping
Tightly coupled to the parameters of Participant()
Definition XMLTag.hpp:24
T visit(T... args)