MappingConfiguration.cpp

Go to the documentation of this file.

#include "MappingConfiguration.hpp"
#include <Eigen/Core>
#include <algorithm>
#include <cstring>
#include <list>
#include <memory>
#include <ostream>
#include <utility>
#include <variant>
#include "logging/LogMacros.hpp"
#include "mapping/AxialGeoMultiscaleMapping.hpp"
#include "mapping/GinkgoRadialBasisFctSolver.hpp"
#include "mapping/LinearCellInterpolationMapping.hpp"
#include "mapping/Mapping.hpp"
#include "mapping/NearestNeighborGradientMapping.hpp"
#include "mapping/NearestNeighborMapping.hpp"
#include "mapping/NearestProjectionMapping.hpp"
#include "mapping/PartitionOfUnityMapping.hpp"
#include "mapping/PetRadialBasisFctMapping.hpp"
#include "mapping/RadialBasisFctMapping.hpp"
#include "mapping/RadialBasisFctSolver.hpp"
#include "mapping/RadialGeoMultiscaleMapping.hpp"
#include "mapping/impl/BasisFunctions.hpp"
#include "mesh/Mesh.hpp"
#include "mesh/SharedPointer.hpp"
#include "mesh/config/MeshConfiguration.hpp"
#include "utils/Parallel.hpp"
#include "utils/Petsc.hpp"
#include "utils/assertion.hpp"
#include "xml/ConfigParser.hpp"
#include "xml/XMLAttribute.hpp"
#include "xml/XMLTag.hpp"
namespace precice::mapping {
 
namespace {
 
// given a list of subtags and parent tags, this function adds all subtags to all
// parent tags
void addSubtagsToParents(std::list<xml::XMLTag> &subtags,
                         std::list<xml::XMLTag> &parents)
{
  for (auto &p : parents) {
    p.addSubtags(subtags);
  }
}
 
// this function uses std::variant in order to add attributes of any type (double, string, bool)
// to all tags in the list of tags \p storage.
using variant_t = std::variant<xml::XMLAttribute<double>, xml::XMLAttribute<std::string>, xml::XMLAttribute<bool>, xml::XMLAttribute<int>>;
template <typename TagStorage>
void addAttributes(TagStorage &storage, const std::vector<variant_t> &attributes)
{
  for (auto &s : storage) {
    for (auto &a : attributes)
      std::visit([&s](auto &&arg) { s.addAttribute(arg); }, a);
  }
}
 
// Enum required for the RBF instantiations
enum struct RBFBackend {
  Eigen,
  PETSc,
  Ginkgo,
  PUM
};
 
// Helper in order to resolve the template instantiations.
// Only the template specializations are of interest
template <RBFBackend Backend, typename RBF>
struct BackendSelector {
  typedef RBF type;
};
 
// Specialization for the RBF Eigen backend
template <typename RBF>
struct BackendSelector<RBFBackend::Eigen, RBF> {
  typedef mapping::RadialBasisFctMapping<RadialBasisFctSolver<RBF>> type;
};
 
// Specialization for the PETSc RBF backend
#ifndef PRECICE_NO_PETSC
template <typename RBF>
struct BackendSelector<RBFBackend::PETSc, RBF> {
  typedef mapping::PetRadialBasisFctMapping<RBF> type;
};
#endif
 
// Specialization for the Ginkgo RBF backend
#ifndef PRECICE_NO_GINKGO
template <typename RBF>
struct BackendSelector<RBFBackend::Ginkgo, RBF> {
  typedef mapping::RadialBasisFctMapping<GinkgoRadialBasisFctSolver<RBF>, MappingConfiguration::GinkgoParameter> type;
};
#endif
// Specialization for the RBF PUM backend
template <typename RBF>
struct BackendSelector<RBFBackend::PUM, RBF> {
  typedef mapping::PartitionOfUnityMapping<RBF> type;
};
 
// Variant holding all available RBF classes
using rbf_variant_t = std::variant<CompactPolynomialC0, CompactPolynomialC2, CompactPolynomialC4, CompactPolynomialC6, CompactPolynomialC8, CompactThinPlateSplinesC2, ThinPlateSplines, VolumeSplines, Multiquadrics, InverseMultiquadrics, Gaussian>;
 
// The actual instantiation of the mapping class, which is called by the visitor \ref getRBFMapping
template <RBFBackend T, typename RADIAL_BASIS_FUNCTION_T, typename... Args>
PtrMapping instantiateRBFMapping(mapping::Mapping::Constraint &constraint, int dimension, RADIAL_BASIS_FUNCTION_T function,
                                 Args &&...args)
{
  return PtrMapping(new typename BackendSelector<T, RADIAL_BASIS_FUNCTION_T>::type(constraint, dimension, function, std::forward<Args>(args)...));
}
 
// Constructs the RBF function based on the functionType
rbf_variant_t constructRBF(BasisFunction functionType, double supportRadius, double shapeParameter)
{
  try {
    switch (functionType) {
    case BasisFunction::WendlandC0: {
      return mapping::CompactPolynomialC0(supportRadius);
    }
    case BasisFunction::WendlandC2: {
      return mapping::CompactPolynomialC2(supportRadius);
    }
    case BasisFunction::WendlandC4: {
      return mapping::CompactPolynomialC4(supportRadius);
    }
    case BasisFunction::WendlandC6: {
      return mapping::CompactPolynomialC6(supportRadius);
    }
    case BasisFunction::WendlandC8: {
      return mapping::CompactPolynomialC8(supportRadius);
    }
    case BasisFunction::CompactThinPlateSplinesC2: {
      return mapping::CompactThinPlateSplinesC2(supportRadius);
    }
    case BasisFunction::ThinPlateSplines: {
      return mapping::ThinPlateSplines();
    }
    case BasisFunction::VolumeSplines: {
      return mapping::VolumeSplines();
    }
    case BasisFunction::Multiquadrics: {
      return mapping::Multiquadrics(shapeParameter);
    }
    case BasisFunction::InverseMultiquadrics: {
      return mapping::InverseMultiquadrics(shapeParameter);
    }
    case BasisFunction::Gaussian: {
      return mapping::Gaussian(shapeParameter);
    }
    default:
      PRECICE_UNREACHABLE("No instantiation was found for the selected basis function.");
    }
  } catch (std::invalid_argument &e) {
    logging::Logger _log{"MappingConfiguration"};
    PRECICE_ERROR(e.what());
  }
}
 
// The actual instantion helper, which avoids enumerating all mapping implementations (more will come) with all RBF kernels
// The first three arguments of the constructor are prescribed: constraint, dimension and the RBF function object, all other
// constructor arguments are just forwarded. The first argument (BasisFunction) indicates then the actual instantiation to return.
template <RBFBackend T, typename... Args>
PtrMapping getRBFMapping(BasisFunction functionType, mapping::Mapping::Constraint &constraint, int dimension, double supportRadius, double shapeParameter,
                         Args &&...args)
{
  // First, construct the RBF function
  auto functionVariant = constructRBF(functionType, supportRadius, shapeParameter);
  // ... and instantiate the corresponding RBF mapping class
  return std::visit([&](auto &&func) { return instantiateRBFMapping<T>(constraint, dimension, func, std::forward<Args>(args)...); }, functionVariant);
}
} // namespace
 
MappingConfiguration::MappingConfiguration(
    xml::XMLTag               &parent,
    mesh::PtrMeshConfiguration meshConfiguration)
    : _meshConfig(std::move(meshConfiguration))
{
  PRECICE_ASSERT(_meshConfig);
  using namespace xml;
 
  // First, we create the available tags
  XMLTag::Occurrence occ = XMLTag::OCCUR_ARBITRARY;
  std::list<XMLTag>  projectionTags{
      XMLTag{*this, TYPE_NEAREST_NEIGHBOR, occ, TAG}.setDocumentation("Nearest-neighbour mapping which uses a rstar-spatial index tree to index meshes and run nearest-neighbour queries."),
      XMLTag{*this, TYPE_NEAREST_PROJECTION, occ, TAG}.setDocumentation("Nearest-projection mapping which uses a rstar-spatial index tree to index meshes and locate the nearest projections."),
      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."),
      XMLTag{*this, TYPE_LINEAR_CELL_INTERPOLATION, occ, TAG}.setDocumentation("Linear cell interpolation mapping which uses a rstar-spatial index tree to index meshes and locate the nearest cell. Only supports 2D meshes.")};
  std::list<XMLTag> rbfDirectTags{
      XMLTag{*this, TYPE_RBF_GLOBAL_DIRECT, occ, TAG}.setDocumentation("Radial-basis-function mapping using a direct solver with a gather-scatter parallelism.")};
  std::list<XMLTag> rbfIterativeTags{
      XMLTag{*this, TYPE_RBF_GLOBAL_ITERATIVE, occ, TAG}.setDocumentation("Radial-basis-function mapping using an iterative solver with a distributed parallelism.")};
  std::list<XMLTag> pumDirectTags{
      XMLTag{*this, TYPE_RBF_PUM_DIRECT, occ, TAG}.setDocumentation("Radial-basis-function mapping using a partition of unity method, which supports a distributed parallelism.")};
  std::list<XMLTag> rbfAliasTag{
      XMLTag{*this, TYPE_RBF_ALIAS, occ, TAG}.setDocumentation("Alias tag, which auto-selects a radial-basis-function mapping depending on the simulation parameter,")};
  std::list<XMLTag> coarseGrainingTags{
      XMLTag{*this, TYPE_COARSE_GRAINING, occ, TAG}.setDocumentation("Coarse graining specifically designed for particle-mesh coupling to write data from the particles to the mesh. The mapping transforms an extensive quantity (e.g., volume, force) into an intensive quantity (e.g., porosity, force-density). "
                                                                     " Currently implemented as just-in-time mapping. Although the constraint does not really fit here (the input is conservative, the output not), we classify it as \"conservative\" for the configuration.")};
  std::list<XMLTag> axialGeoMultiscaleTags{
      XMLTag{*this, TYPE_AXIAL_GEOMETRIC_MULTISCALE, occ, TAG}
          .setDocumentation("Axial geometric multiscale mapping between one 1D and multiple 3D vertices.")};
  std::list<XMLTag> radialGeoMultiscaleTags{
      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.")};
 
  // List of all attributes with corresponding documentation
  auto attrDirection = XMLAttribute<std::string>(ATTR_DIRECTION)
                           .setOptions({DIRECTION_WRITE, DIRECTION_READ})
                           .setDocumentation("Write mappings map written data prior to communication, thus in the same participant who writes the data. "
                                             "Read mappings map received data after communication, thus in the same participant who reads the data.");
 
  auto attrFromMesh = XMLAttribute<std::string>(ATTR_FROM, "")
                          .setDocumentation("The mesh to map the data from. The default name is an empty mesh name, which is only valid for a just-in-time mapping (using the API functions \"writeAndMapData\" or \"mapAndReadData\").");
 
  auto attrToMesh = XMLAttribute<std::string>(ATTR_TO, "")
                        .setDocumentation("The mesh to map the data to. The default name is an empty mesh name, which is only valid for a just-in-time mapping (using the API functions \"writeAndMapData\" or \"mapAndReadData\").");
 
  auto attrConstraint = XMLAttribute<std::string>(ATTR_CONSTRAINT)
                            .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.")
                            .setOptions({CONSTRAINT_CONSERVATIVE, CONSTRAINT_CONSISTENT, CONSTRAINT_SCALED_CONSISTENT_SURFACE, CONSTRAINT_SCALED_CONSISTENT_VOLUME});
  auto attrXDead = makeXMLAttribute(ATTR_X_DEAD, false)
                       .setDocumentation("If set to true, the x axis will be ignored for the mapping");
  auto attrYDead = makeXMLAttribute(ATTR_Y_DEAD, false)
                       .setDocumentation("If set to true, the y axis will be ignored for the mapping");
  auto attrZDead = makeXMLAttribute(ATTR_Z_DEAD, false)
                       .setDocumentation("If set to true, the z axis will be ignored for the mapping");
  auto attrPolynomial = makeXMLAttribute(ATTR_POLYNOMIAL, POLYNOMIAL_SEPARATE)
                            .setDocumentation("Toggles use of the global polynomial")
                            .setOptions({POLYNOMIAL_ON, POLYNOMIAL_OFF, POLYNOMIAL_SEPARATE});
  auto attrPumPolynomial = makeXMLAttribute(ATTR_POLYNOMIAL, POLYNOMIAL_SEPARATE)
                               .setDocumentation("Toggles use a local (per cluster) polynomial")
                               .setOptions({POLYNOMIAL_OFF, POLYNOMIAL_SEPARATE});
 
  auto attrcgRadius = makeXMLAttribute<double>(ATTR_CG_RADIUS, 0.)
                          .setDocumentation("Radius or range of the coarsening function (Lucy function).");
 
  auto attrSolverRtol = makeXMLAttribute(ATTR_SOLVER_RTOL, 1e-9)
                            .setDocumentation("Solver relative tolerance for convergence");
  // TODO: Discuss whether we wanto to introduce this attribute
  // auto attrMaxIterations = makeXMLAttribute(ATTR_MAX_ITERATIONS, 1e6)
  //                              .setDocumentation("Maximum number of iterations of the solver");
 
  auto verticesPerCluster = XMLAttribute<int>(ATTR_VERTICES_PER_CLUSTER, 50)
                                .setDocumentation("Average number of vertices per cluster (partition) applied in the rbf partition of unity method.");
  auto relativeOverlap = makeXMLAttribute(ATTR_RELATIVE_OVERLAP, 0.15)
                             .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.");
  auto projectToInput = XMLAttribute<bool>(ATTR_PROJECT_TO_INPUT, true)
                            .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.");
 
  auto attrGeoMultiscaleDimension = XMLAttribute<std::string>(ATTR_GEOMETRIC_MULTISCALE_DIMENSION)
                                        .setDocumentation("Specifies the dimensionality pairing used in geometric multiscale mapping. Options: '1D-3D', '1D-2D' or '2D-3D'.")
                                        .setOptions({GEOMETRIC_MULTISCALE_DIMENSION_1D3D, GEOMETRIC_MULTISCALE_DIMENSION_1D2D, GEOMETRIC_MULTISCALE_DIMENSION_2D3D});
  auto attrGeoMultiscaleType = XMLAttribute<std::string>(ATTR_GEOMETRIC_MULTISCALE_TYPE)
                                   .setDocumentation("Type of geometric multiscale mapping. Either 'spread' or 'collect'.")
                                   .setOptions({GEOMETRIC_MULTISCALE_TYPE_SPREAD, GEOMETRIC_MULTISCALE_TYPE_COLLECT});
  auto attrGeoMultiscaleAxis = XMLAttribute<std::string>(ATTR_GEOMETRIC_MULTISCALE_AXIS)
                                   .setDocumentation("Principle axis along which geometric multiscale mapping is performed.")
                                   .setOptions({GEOMETRIC_MULTISCALE_AXIS_X, GEOMETRIC_MULTISCALE_AXIS_Y, GEOMETRIC_MULTISCALE_AXIS_Z});
  auto attrGeoMultiscaleRadius = XMLAttribute<double>(ATTR_GEOMETRIC_MULTISCALE_RADIUS)
                                     .setDocumentation("Radius of the cross-sectional interface between the participants.");
  auto attrGeoMultiscaleCrossSectionProfile = XMLAttribute<std::string>(ATTR_GEOMETRIC_MULTISCALE_CROSS_SECTION_PROFILE)
                                                  .setDocumentation("Profile of the mapped variable along the cross-sectional interface: 'uniform' or 'parabolic'")
                                                  .setOptions({GEOMETRIC_MULTISCALE_CROSS_SECTION_PROFILE_UNIFORM, GEOMETRIC_MULTISCALE_CROSS_SECTION_PROFILE_PARABOLIC})
                                                  .setDefaultValue(GEOMETRIC_MULTISCALE_CROSS_SECTION_PROFILE_UNIFORM);
  auto attrGeoMultiscaleCrossSection = XMLAttribute<std::string>(ATTR_GEOMETRIC_MULTISCALE_CROSS_SECTION)
                                           .setDocumentation("Cross section of the interface of the participants: 'circle' or 'square'")
                                           .setOptions({GEOMETRIC_MULTISCALE_CROSS_SECTION_CIRCLE, GEOMETRIC_MULTISCALE_CROSS_SECTION_SQUARE})
                                           .setDefaultValue(GEOMETRIC_MULTISCALE_CROSS_SECTION_CIRCLE);
 
  // Add the relevant attributes to the relevant tags
  addAttributes(projectionTags, {attrFromMesh, attrToMesh, attrDirection, attrConstraint});
  addAttributes(rbfDirectTags, {attrFromMesh, attrToMesh, attrDirection, attrConstraint, attrPolynomial, attrXDead, attrYDead, attrZDead});
  addAttributes(rbfIterativeTags, {attrFromMesh, attrToMesh, attrDirection, attrConstraint, attrPolynomial, attrXDead, attrYDead, attrZDead, attrSolverRtol});
  addAttributes(pumDirectTags, {attrFromMesh, attrToMesh, attrDirection, attrConstraint, attrPumPolynomial, verticesPerCluster, relativeOverlap, projectToInput});
  addAttributes(rbfAliasTag, {attrFromMesh, attrToMesh, attrDirection, attrConstraint, attrXDead, attrYDead, attrZDead});
  addAttributes(coarseGrainingTags, {attrFromMesh, attrToMesh, attrDirection, attrConstraint, attrcgRadius});
  addAttributes(axialGeoMultiscaleTags, {attrFromMesh, attrToMesh, attrDirection, attrConstraint});
  addAttributes(radialGeoMultiscaleTags, {attrFromMesh, attrToMesh, attrDirection, attrConstraint});
  addAttributes(axialGeoMultiscaleTags, {attrGeoMultiscaleDimension, attrGeoMultiscaleType, attrGeoMultiscaleAxis, attrGeoMultiscaleRadius, attrGeoMultiscaleCrossSectionProfile, attrGeoMultiscaleCrossSection});
  addAttributes(radialGeoMultiscaleTags, {attrGeoMultiscaleType, attrGeoMultiscaleAxis, attrGeoMultiscaleRadius});
 
  // Now we take care of the subtag executor. We repeat some of the subtags in order to add individual documentation
  XMLTag::Occurrence once = XMLTag::OCCUR_NOT_OR_ONCE;
  // TODO, make type an int
  auto attrDeviceId = makeXMLAttribute(ATTR_DEVICE_ID, static_cast<std::string>("0"))
                          .setDocumentation("Setting of the GPU device: Set \"auto\" to assign GPUs to each MPI rank in a round robin fashion or specify a number between 0 and the number of available GPUs-1 to assign all MPI ranks to one GPU device with the given ID.");
  auto attrNThreads = makeXMLAttribute(ATTR_N_THREADS, static_cast<int>(0))
                          .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.");
 
  // First, we have the executors for the direct solvers
  {
    std::list<XMLTag> cpuExecutor{
        XMLTag{*this, EXECUTOR_CPU, once, SUBTAG_EXECUTOR}.setDocumentation("The default executor, which uses a single-core CPU with a gather-scatter parallelism.")};
    std::list<XMLTag> deviceExecutors{
        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."),
        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.")};
 
    addAttributes(deviceExecutors, {attrDeviceId});
    addSubtagsToParents(cpuExecutor, rbfDirectTags);
    addSubtagsToParents(deviceExecutors, rbfDirectTags);
  }
  // Second, the executors for the iterative solver
  {
    std::list<XMLTag> cpuExecutor{
        XMLTag{*this, EXECUTOR_CPU, once, SUBTAG_EXECUTOR}.setDocumentation("The default executor relying on PETSc, which uses CPUs and distributed memory parallelism via MPI.")};
    std::list<XMLTag> deviceExecutors{
        XMLTag{*this, EXECUTOR_CUDA, once, SUBTAG_EXECUTOR}.setDocumentation("Cuda (Nvidia) executor, which uses Ginkgo with a gather-scatter parallelism."),
        XMLTag{*this, EXECUTOR_HIP, once, SUBTAG_EXECUTOR}.setDocumentation("Hip (AMD/Nvidia) executor, which uses hipSolver with a gather-scatter parallelism.")};
    std::list<XMLTag> ompExecutor{
        XMLTag{*this, EXECUTOR_OMP, once, SUBTAG_EXECUTOR}.setDocumentation("OpenMP executor, which uses Ginkgo with a gather-scatter parallelism.")};
 
    addAttributes(deviceExecutors, {attrDeviceId});
    addAttributes(ompExecutor, {attrNThreads});
    addSubtagsToParents(cpuExecutor, rbfIterativeTags);
    addSubtagsToParents(deviceExecutors, rbfIterativeTags);
    addSubtagsToParents(ompExecutor, rbfIterativeTags);
  }
  {
    std::list<XMLTag> cpuExecutor{
        XMLTag{*this, EXECUTOR_CPU, once, SUBTAG_EXECUTOR}.setDocumentation("The default executor using a CPU and a distributed memory parallelism via MPI.")};
    std::list<XMLTag> deviceExecutors{
        XMLTag{*this, EXECUTOR_CUDA, once, SUBTAG_EXECUTOR}.setDocumentation("Cuda (Nvidia) executor, which uses Kokkos-kernels, fully parallel"),
        XMLTag{*this, EXECUTOR_HIP, once, SUBTAG_EXECUTOR}.setDocumentation("Hip (AMD/Nvidia) executor, which uses Kokkos-kernels, fully parallel."),
        XMLTag{*this, EXECUTOR_SYCL, once, SUBTAG_EXECUTOR}.setDocumentation("SYCL (e.g. Intel) executor, which uses Kokkos-kernels, fully parallel.")};
    std::list<XMLTag> ompExecutor{
        XMLTag{*this, EXECUTOR_OMP, once, SUBTAG_EXECUTOR}.setDocumentation("OpenMP executor, which uses Kokkos-kernel, fully parallel.")};
 
    auto attrExecutionMode = makeXMLAttribute(ATTR_EXECUTION_MODE, "minimal-memory")
                                 .setDocumentation("Toggle to switch between a minimal-memory vs a minimal-compute algorithm. For option \"minimal-memory\", the RBF evaluation is recomputed for each data mapping on-the-fly, which saves approximately half the memory consumption (if meshes have a similar resolution), but may (!) be slower (depends heavily on the hardware). "
                                                   "For the option \"minimal-compute\", the RBF evaluation is precomputed, which may be faster, but consumes more memory.")
                                 .setOptions({"minimal-memory", "minimal-compute"});
 
    addAttributes(deviceExecutors, {attrDeviceId, attrExecutionMode});
    addAttributes(ompExecutor, {attrNThreads, attrExecutionMode});
    addSubtagsToParents(cpuExecutor, pumDirectTags);
    addSubtagsToParents(deviceExecutors, pumDirectTags);
    addSubtagsToParents(ompExecutor, pumDirectTags);
  }
  // The alias tag doesn't receive the subtag at all
 
  // Now we take care of the subtag basis function
  // First, we have the tags using a support radius
  std::list<XMLTag> supportRadiusRBF{
      XMLTag{*this, RBF_CPOLYNOMIAL_C0, once, SUBTAG_BASIS_FUNCTION}.setDocumentation("Wendland C0 function"),
      XMLTag{*this, RBF_CPOLYNOMIAL_C2, once, SUBTAG_BASIS_FUNCTION}.setDocumentation("Wendland C2 function"),
      XMLTag{*this, RBF_CPOLYNOMIAL_C4, once, SUBTAG_BASIS_FUNCTION}.setDocumentation("Wendland C4 function"),
      XMLTag{*this, RBF_CPOLYNOMIAL_C6, once, SUBTAG_BASIS_FUNCTION}.setDocumentation("Wendland C6 function"),
      XMLTag{*this, RBF_CPOLYNOMIAL_C8, once, SUBTAG_BASIS_FUNCTION}.setDocumentation("Wendland C8 function"),
      XMLTag{*this, RBF_CTPS_C2, once, SUBTAG_BASIS_FUNCTION}.setDocumentation("Compact thin-plate-spline C2")};
 
  auto attrSupportRadius = XMLAttribute<double>(ATTR_SUPPORT_RADIUS)
                               .setDocumentation("Support radius of each RBF basis function (global choice).");
 
  addAttributes(supportRadiusRBF, {attrSupportRadius});
  addSubtagsToParents(supportRadiusRBF, rbfIterativeTags);
  addSubtagsToParents(supportRadiusRBF, rbfDirectTags);
  addSubtagsToParents(supportRadiusRBF, pumDirectTags);
  addSubtagsToParents(supportRadiusRBF, rbfAliasTag);
 
  // Now the tags using a shape parameter
  std::list<XMLTag> shapeParameterRBF{
      XMLTag{*this, RBF_MULTIQUADRICS, once, SUBTAG_BASIS_FUNCTION}.setDocumentation("Multiquadrics"),
      XMLTag{*this, RBF_INV_MULTIQUADRICS, once, SUBTAG_BASIS_FUNCTION}.setDocumentation("Inverse multiquadrics")};
 
  auto attrShapeParam = XMLAttribute<double>(ATTR_SHAPE_PARAM)
                            .setDocumentation("Specific shape parameter for RBF basis function.");
 
  addAttributes(shapeParameterRBF, {attrShapeParam});
  addSubtagsToParents(shapeParameterRBF, rbfIterativeTags);
  addSubtagsToParents(shapeParameterRBF, rbfDirectTags);
  addSubtagsToParents(shapeParameterRBF, pumDirectTags);
  addSubtagsToParents(shapeParameterRBF, rbfAliasTag);
 
  // For the Gaussian, we need default values as the user can pass a support radius or a shape parameter
  std::list<XMLTag> GaussRBF{
      XMLTag{*this, RBF_GAUSSIAN, once, SUBTAG_BASIS_FUNCTION}.setDocumentation("Gaussian basis function accepting a support radius or a shape parameter.")};
  attrShapeParam.setDefaultValue(std::numeric_limits<double>::quiet_NaN());
  attrSupportRadius.setDefaultValue(std::numeric_limits<double>::quiet_NaN());
  addAttributes(GaussRBF, {attrShapeParam, attrSupportRadius});
  addSubtagsToParents(GaussRBF, rbfIterativeTags);
  addSubtagsToParents(GaussRBF, rbfDirectTags);
  addSubtagsToParents(GaussRBF, pumDirectTags);
  addSubtagsToParents(GaussRBF, rbfAliasTag);
 
  // tags without an attribute
  std::list<XMLTag> attributelessRBFs{
      XMLTag{*this, RBF_TPS, once, SUBTAG_BASIS_FUNCTION}.setDocumentation("Thin-plate-splines"),
      XMLTag{*this, RBF_VOLUME_SPLINES, once, SUBTAG_BASIS_FUNCTION}.setDocumentation("Volume splines")};
 
  addSubtagsToParents(attributelessRBFs, rbfIterativeTags);
  addSubtagsToParents(attributelessRBFs, rbfDirectTags);
  addSubtagsToParents(attributelessRBFs, pumDirectTags);
  addSubtagsToParents(attributelessRBFs, rbfAliasTag);
 
  // Add all tags to the mapping tag
  parent.addSubtags(projectionTags);
  parent.addSubtags(rbfIterativeTags);
  parent.addSubtags(rbfDirectTags);
  parent.addSubtags(pumDirectTags);
  parent.addSubtags(rbfAliasTag);
  parent.addSubtags(coarseGrainingTags);
  parent.addSubtags(axialGeoMultiscaleTags);
  parent.addSubtags(radialGeoMultiscaleTags);
}
 
void MappingConfiguration::setExperimental(
    bool experimental)
{
  _experimental = experimental;
}
 
void MappingConfiguration::xmlTagCallback(
    const xml::ConfigurationContext &context,
    xml::XMLTag                     &tag)
{
  PRECICE_TRACE(tag.getName());
  if (tag.getNamespace() == TAG) {
    // Mandatory tags
    std::string dir        = tag.getStringAttributeValue(ATTR_DIRECTION);
    std::string fromMesh   = tag.getStringAttributeValue(ATTR_FROM);
    std::string toMesh     = tag.getStringAttributeValue(ATTR_TO);
    std::string type       = tag.getName();
    std::string constraint = tag.getStringAttributeValue(ATTR_CONSTRAINT);
 
    // Check that either of the two is provided
    PRECICE_CHECK(!toMesh.empty() || !fromMesh.empty(), "Neither a \"to\" nor a \"from\" mesh was defined in a preCICE mapping. Most data mappings require a \"to\" and a \"from\" mesh. "
                                                        "For just-in-time mapping, at least one of both attributes has to be specified.");
 
    // Restrict to read-consistent and write-conservative for just-in-time mapping
 
    // The from mesh cannot be empty due to the check above
    PRECICE_CHECK(!toMesh.empty() || (toMesh.empty() && dir == DIRECTION_READ && constraint == CONSTRAINT_CONSISTENT),
                  "The mapping from mesh \"{0}\" has no \"to\" mesh, which configures a just-in-time mapping for the mesh \"{0}\". For just-in-time mapping, only read-consistent (direction = \"read\" and no \"to\" mesh) and write-conservative (direction = \"write\" and no \"from\" mesh) are implemented.", fromMesh);
    PRECICE_CHECK(!fromMesh.empty() || (fromMesh.empty() && dir == DIRECTION_WRITE && constraint == CONSTRAINT_CONSERVATIVE),
                  "The mapping to mesh \"{0}\" has no \"from\" mesh, which configures a just-in-time mapping for the mesh \"{0}\". For just-in-time mapping, only read-consistent (direction = \"read\" and no \"to\" mesh) and write-conservative (direction = \"write\" and no \"from\" mesh) are implemented.", toMesh);
 
    PRECICE_INFO_IF(toMesh.empty(), "Using just-in-time mapping from mesh \"{}\"", fromMesh);
    PRECICE_INFO_IF(fromMesh.empty(), "Using just-in-time mapping to mesh \"{}\"", toMesh);
 
    // optional tags
    // We set here default values, but their actual value doesn't really matter.
    // It's just for the mapping methods, which do not use these attributes at all.
    bool        xDead         = tag.getBooleanAttributeValue(ATTR_X_DEAD, false);
    bool        yDead         = tag.getBooleanAttributeValue(ATTR_Y_DEAD, false);
    bool        zDead         = tag.getBooleanAttributeValue(ATTR_Z_DEAD, false);
    double      cgRadius      = tag.getDoubleAttributeValue(ATTR_CG_RADIUS, 0.);
    double      solverRtol    = tag.getDoubleAttributeValue(ATTR_SOLVER_RTOL, 1e-9);
    std::string strPolynomial = tag.getStringAttributeValue(ATTR_POLYNOMIAL, POLYNOMIAL_SEPARATE);
 
    // geometric multiscale related tags
    std::string geoMultiscaleDimension    = tag.getStringAttributeValue(ATTR_GEOMETRIC_MULTISCALE_DIMENSION, "");
    std::string geoMultiscaleType         = tag.getStringAttributeValue(ATTR_GEOMETRIC_MULTISCALE_TYPE, "");
    std::string geoMultiscaleAxis         = tag.getStringAttributeValue(ATTR_GEOMETRIC_MULTISCALE_AXIS, "");
    double      multiscaleRadius          = tag.getDoubleAttributeValue(ATTR_GEOMETRIC_MULTISCALE_RADIUS, 1.0);
    std::string geoMultiscaleProfile      = tag.getStringAttributeValue(ATTR_GEOMETRIC_MULTISCALE_CROSS_SECTION_PROFILE, "");
    std::string geoMultiscaleCrossSection = tag.getStringAttributeValue(ATTR_GEOMETRIC_MULTISCALE_CROSS_SECTION, "");
 
    if (type == TYPE_AXIAL_GEOMETRIC_MULTISCALE || type == TYPE_RADIAL_GEOMETRIC_MULTISCALE || type == TYPE_COARSE_GRAINING) {
      PRECICE_CHECK(_experimental, "The configured mapping \"{}\" is experimental and the configuration can change between minor releases. Set experimental=\"on\" in the precice-configuration tag.", type);
    }
 
    if (type == TYPE_AXIAL_GEOMETRIC_MULTISCALE && context.size > 1) {
      throw std::runtime_error{"Axial geometric multiscale mapping is not available for parallel participants."};
    }
 
    if (type == TYPE_RADIAL_GEOMETRIC_MULTISCALE && context.size > 1) {
      throw std::runtime_error{"Radial geometric multiscale mapping is not available for parallel participants."};
    }
 
    // pum related tags
    int    verticesPerCluster = tag.getIntAttributeValue(ATTR_VERTICES_PER_CLUSTER, 100);
    double relativeOverlap    = tag.getDoubleAttributeValue(ATTR_RELATIVE_OVERLAP, 0.3);
    bool   projectToInput     = tag.getBooleanAttributeValue(ATTR_PROJECT_TO_INPUT, true);
 
    // Convert raw string into enum types as the constructors take enums
    if (constraint == CONSTRAINT_CONSERVATIVE) {
      constraintValue = Mapping::CONSERVATIVE;
    } else if (constraint == CONSTRAINT_CONSISTENT) {
      constraintValue = Mapping::CONSISTENT;
    } else if (constraint == CONSTRAINT_SCALED_CONSISTENT_SURFACE) {
      constraintValue = Mapping::SCALED_CONSISTENT_SURFACE;
    } else if (constraint == CONSTRAINT_SCALED_CONSISTENT_VOLUME) {
      constraintValue = Mapping::SCALED_CONSISTENT_VOLUME;
    } else {
      PRECICE_UNREACHABLE("Unknown mapping constraint \"{}\".", constraint);
    }
 
    ConfiguredMapping configuredMapping = createMapping(dir, type, fromMesh, toMesh, cgRadius, geoMultiscaleDimension, geoMultiscaleType, geoMultiscaleAxis, multiscaleRadius, geoMultiscaleProfile, geoMultiscaleCrossSection);
 
    _rbfConfig = configureRBFMapping(type, strPolynomial, xDead, yDead, zDead, solverRtol, verticesPerCluster, relativeOverlap, projectToInput);
 
    checkDuplicates(configuredMapping);
    _mappings.push_back(configuredMapping);
  } else if (tag.getNamespace() == SUBTAG_BASIS_FUNCTION) {
 
    PRECICE_ASSERT(!_mappings.empty());
    PRECICE_CHECK(_mappings.back().requiresBasisFunction == true, "A basis-function was defined for the mapping "
                                                                  "from mesh \"{}\" to mesh \"{}\", but no basis-function is required for this mapping type. "
                                                                  "Please remove the basis-function tag or configure an rbf mapping, which requires a basis-function.",
                  _mappings.back().fromMesh->getName(), _mappings.back().toMesh->getName());
    // We can only set one subtag
    PRECICE_CHECK(_rbfConfig.basisFunctionDefined == false, "More than one basis-function was defined for the mapping "
                                                            "from mesh \"{}\" to mesh \"{}\".",
                  _mappings.back().fromMesh->getName(), _mappings.back().toMesh->getName());
 
    std::string basisFctName   = tag.getName();
    double      supportRadius  = tag.getDoubleAttributeValue(ATTR_SUPPORT_RADIUS, 0.);
    double      shapeParameter = tag.getDoubleAttributeValue(ATTR_SHAPE_PARAM, 0.);
 
    _rbfConfig.basisFunction        = parseBasisFunctions(basisFctName);
    _rbfConfig.basisFunctionDefined = true;
    // The Gaussian RBF is always treated as a shape-parameter RBF. Hence, we have to convert the support radius, if necessary
    if (_rbfConfig.basisFunction == BasisFunction::Gaussian) {
      const bool exactlyOneSet = (std::isfinite(supportRadius) && !std::isfinite(shapeParameter)) ||
                                 (std::isfinite(shapeParameter) && !std::isfinite(supportRadius));
      PRECICE_CHECK(exactlyOneSet, "The specified parameters for the Gaussian RBF mapping are invalid. Please specify either a \"shape-parameter\" or a \"support-radius\".");
 
      if (std::isfinite(supportRadius) && !std::isfinite(shapeParameter)) {
        shapeParameter = std::sqrt(-std::log(Gaussian::cutoffThreshold)) / supportRadius;
      }
    }
 
    _rbfConfig.supportRadius  = supportRadius;
    _rbfConfig.shapeParameter = shapeParameter;
  } else if (tag.getNamespace() == SUBTAG_EXECUTOR) {
    _executorConfig = std::make_unique<ExecutorConfiguration>();
 
    if (tag.getName() == EXECUTOR_CPU) {
      _executorConfig->executor = ExecutorConfiguration::Executor::CPU;
    } else if (tag.getName() == EXECUTOR_CUDA) {
      _executorConfig->executor = ExecutorConfiguration::Executor::CUDA;
    } else if (tag.getName() == EXECUTOR_HIP) {
      _executorConfig->executor = ExecutorConfiguration::Executor::HIP;
    } else if (tag.getName() == EXECUTOR_SYCL) {
      _executorConfig->executor = ExecutorConfiguration::Executor::SYCL;
    } else if (tag.getName() == EXECUTOR_OMP) {
      _executorConfig->executor = ExecutorConfiguration::Executor::OpenMP;
    }
 
    auto did = tag.getStringAttributeValue(ATTR_DEVICE_ID, "0");
    if (did == "auto") {
      _executorConfig->deviceId = -1;
    } else {
      try {
        _executorConfig->deviceId = std::stoi(did);
        PRECICE_CHECK(_executorConfig->deviceId >= 0, "The argument provided to \"gpu-device-id\" in the precice configuration file is invalid (negative device id)");
      } catch (const std::invalid_argument &e) {
        throw precice::Error("The argument provided to \"gpu-device-id\" in the precice configuration file is invalid (not a valid input).");
      } catch (const std::out_of_range &e) {
        throw precice::Error("The argument provided to \"gpu-device-id\" in the precice configuration file is invalid (out of range).");
      }
    }
 
    _executorConfig->nThreads = tag.getIntAttributeValue(ATTR_N_THREADS, 0);
    auto mode                 = tag.getStringAttributeValue(ATTR_EXECUTION_MODE, "minimal-compute");
    if (mode == "minimal-memory")
      _executorConfig->computeEvaluationOffline = false;
    else if (mode == "minimal-compute")
      _executorConfig->computeEvaluationOffline = true;
    else {
      PRECICE_UNREACHABLE("Unknown execution mode");
    }
  }
}
 
MappingConfiguration::RBFConfiguration MappingConfiguration::configureRBFMapping(const std::string &type,
                                                                                 const std::string &polynomial,
                                                                                 bool xDead, bool yDead, bool zDead,
                                                                                 double solverRtol,
                                                                                 double verticesPerCluster,
                                                                                 double relativeOverlap,
                                                                                 bool   projectToInput) const
{
  RBFConfiguration rbfConfig;
 
  if (type == TYPE_RBF_GLOBAL_ITERATIVE)
    rbfConfig.solver = RBFConfiguration::SystemSolver::GlobalIterative;
  else if (type == TYPE_RBF_GLOBAL_DIRECT)
    rbfConfig.solver = RBFConfiguration::SystemSolver::GlobalDirect;
  else if (type == TYPE_RBF_PUM_DIRECT)
    rbfConfig.solver = RBFConfiguration::SystemSolver::PUMDirect;
  else {
    // Rather simple auto-selection (for now), consisting of the PUM Eigen backend
    rbfConfig.solver = RBFConfiguration::SystemSolver::PUMDirect;
 
    // A more sophisticated criterion here could take the globalNumberOfVertices into account
    // (the mesh pointer is stored in the configuredMapping anyway), but this quantity is not yet computed
    // during the configuration time.
  }
 
  if (polynomial == POLYNOMIAL_SEPARATE)
    rbfConfig.polynomial = Polynomial::SEPARATE;
  else if (polynomial == POLYNOMIAL_ON)
    rbfConfig.polynomial = Polynomial::ON;
  else if (polynomial == POLYNOMIAL_OFF)
    rbfConfig.polynomial = Polynomial::OFF;
  else
    PRECICE_UNREACHABLE("Unknown polynomial configuration.");
 
  rbfConfig.deadAxis   = {{xDead, yDead, zDead}};
  rbfConfig.solverRtol = solverRtol;
 
  rbfConfig.verticesPerCluster = verticesPerCluster;
  rbfConfig.relativeOverlap    = relativeOverlap;
  rbfConfig.projectToInput     = projectToInput;
 
  return rbfConfig;
}
 
MappingConfiguration::ConfiguredMapping MappingConfiguration::createMapping(
    const std::string &direction,
    const std::string &type,
    const std::string &fromMeshName,
    const std::string &toMeshName,
    const double       cgRadius,
    const std::string &geoMultiscaleDimension,
    const std::string &geoMultiscaleType,
    const std::string &geoMultiscaleAxis,
    const double      &multiscaleRadius,
    const std::string &geoMultiscaleProfile,
    const std::string &geoMultiscaleCrossSection) const
{
  PRECICE_TRACE(direction, type);
 
  ConfiguredMapping configuredMapping;
  mesh::PtrMesh     fromMesh(_meshConfig->getMesh(fromMeshName));
  mesh::PtrMesh     toMesh(_meshConfig->getMesh(toMeshName));
 
  // Handle for just-in-time mapping, we copy over the dimension and leave everything else
  if (!fromMesh && fromMeshName.empty()) {
    PRECICE_CHECK(toMesh,
                  "Mesh \"{0}\" was not found while creating a mapping. "
                  "Please correct the to=\"{0}\" attribute.",
                  toMeshName);
    fromMesh = mesh::MeshConfiguration::getJustInTimeMappingMesh(toMesh->getDimensions());
  }
  if (!toMesh && toMeshName.empty()) {
    PRECICE_CHECK(fromMesh,
                  "Mesh \"{0}\" was not found while creating a mapping. "
                  "Please correct the from=\"{0}\" attribute.",
                  fromMeshName);
    toMesh = mesh::MeshConfiguration::getJustInTimeMappingMesh(fromMesh->getDimensions());
  }
 
  PRECICE_CHECK(fromMesh,
                "Mesh \"{0}\" was not found while creating a mapping. "
                "Please correct the from=\"{0}\" attribute.",
                fromMeshName);
  PRECICE_CHECK(toMesh,
                "Mesh \"{0}\" was not found while creating a mapping. "
                "Please correct the to=\"{0}\" attribute.",
                toMeshName);
 
  PRECICE_CHECK((!toMesh->isJustInTime() && !fromMesh->isJustInTime()) ||
                    ((toMesh->isJustInTime() || fromMesh->isJustInTime()) && (type == TYPE_NEAREST_NEIGHBOR || type == TYPE_RBF_PUM_DIRECT || type == TYPE_RBF_ALIAS || type == TYPE_COARSE_GRAINING)),
                "A just-in-time mapping was configured from mesh \"{}\" to mesh \"{}\" using \"mapping:{}\", which is currently not implemented. "
                "Available mapping types are \"mapping:{}\", \"mapping:{}\", \"mapping:{}\" and  \"mapping:{}\".",
                fromMesh->getName(), toMesh->getName(), type, TYPE_NEAREST_NEIGHBOR, TYPE_RBF_ALIAS, TYPE_RBF_PUM_DIRECT, TYPE_COARSE_GRAINING);
 
  // Check for compatible mesh dimensions
  PRECICE_CHECK(fromMesh->getDimensions() == toMesh->getDimensions(),
                "Mapping between meshes of different dimensions is not allowed yet. "
                "Please set the same dimensions attribute to meshes \"{}\" and \"{}\", "
                "or choose different meshes.",
                fromMeshName, toMeshName);
 
  configuredMapping.fromMesh = fromMesh;
  configuredMapping.toMesh   = toMesh;
 
  if (direction == DIRECTION_WRITE) {
    configuredMapping.direction = WRITE;
  } else if (direction == DIRECTION_READ) {
    configuredMapping.direction = READ;
  } else {
    PRECICE_UNREACHABLE("Unknown mapping direction type \"{}\".", direction);
  }
 
  // Create all projection based mappings
  if (type == TYPE_NEAREST_NEIGHBOR) {
    configuredMapping.mapping = PtrMapping(new NearestNeighborMapping(constraintValue, fromMesh->getDimensions()));
  } else if (type == TYPE_NEAREST_PROJECTION) {
    configuredMapping.mapping = PtrMapping(new NearestProjectionMapping(constraintValue, fromMesh->getDimensions()));
  } else if (type == TYPE_LINEAR_CELL_INTERPOLATION) {
    configuredMapping.mapping = PtrMapping(new LinearCellInterpolationMapping(constraintValue, fromMesh->getDimensions()));
  } else if (type == TYPE_COARSE_GRAINING) {
    configuredMapping.mapping = PtrMapping(new CoarseGrainingMapping(constraintValue, fromMesh->getDimensions(), cgRadius));
  } else if (type == TYPE_NEAREST_NEIGHBOR_GRADIENT) {
 
    // NNG is not applicable with the conservative constraint
    PRECICE_CHECK(constraintValue != Mapping::CONSERVATIVE,
                  "Nearest-neighbor-gradient mapping is not implemented using a \"conservative\" constraint. "
                  "Please select constraint=\" consistent\" or a different mapping method.");
 
    configuredMapping.mapping = PtrMapping(new NearestNeighborGradientMapping(constraintValue, fromMesh->getDimensions()));
 
  } else if (type == TYPE_AXIAL_GEOMETRIC_MULTISCALE) {
 
    // Axial geometric multiscale is not applicable with the conservative constraint
    PRECICE_CHECK(constraintValue != Mapping::CONSERVATIVE,
                  "Axial geometric multiscale mapping is not implemented for the \"conservative\" constraint. "
                  "Please select constraint=\" consistent\" or a different mapping method.");
 
    // Convert strings into enums
    AxialGeoMultiscaleMapping::MultiscaleAxis multiscaleAxis;
    if (geoMultiscaleAxis == "x") {
      multiscaleAxis = AxialGeoMultiscaleMapping::MultiscaleAxis::X;
    } else if (geoMultiscaleAxis == "y") {
      multiscaleAxis = AxialGeoMultiscaleMapping::MultiscaleAxis::Y;
    } else if (geoMultiscaleAxis == "z") {
      multiscaleAxis = AxialGeoMultiscaleMapping::MultiscaleAxis::Z;
    } else {
      PRECICE_UNREACHABLE("Unknown geometric multiscale axis \"{}\".", geoMultiscaleAxis);
    }
 
    AxialGeoMultiscaleMapping::MultiscaleDimension multiscaleDimension;
    if (geoMultiscaleDimension == "1d-3d") {
      multiscaleDimension = AxialGeoMultiscaleMapping::MultiscaleDimension::D1D3;
    } else if (geoMultiscaleDimension == "1d-2d") {
      multiscaleDimension = AxialGeoMultiscaleMapping::MultiscaleDimension::D1D2;
    } else if (geoMultiscaleDimension == "2d-3d") {
      multiscaleDimension = AxialGeoMultiscaleMapping::MultiscaleDimension::D2D3;
    } else {
      PRECICE_UNREACHABLE("Unknown dimension \"{}\".", geoMultiscaleDimension);
    }
 
    AxialGeoMultiscaleMapping::MultiscaleType multiscaleType;
    if (geoMultiscaleType == "spread") {
      multiscaleType = AxialGeoMultiscaleMapping::MultiscaleType::SPREAD;
    } else if (geoMultiscaleType == "collect") {
      multiscaleType = AxialGeoMultiscaleMapping::MultiscaleType::COLLECT;
    } else {
      PRECICE_UNREACHABLE("Unknown geometric multiscale type \"{}\".", geoMultiscaleType);
    }
 
    AxialGeoMultiscaleMapping::MultiscaleProfile multiscaleProfile = AxialGeoMultiscaleMapping::MultiscaleProfile::PARABOLIC;
    if (geoMultiscaleProfile == "parabolic") {
      multiscaleProfile = AxialGeoMultiscaleMapping::MultiscaleProfile::PARABOLIC;
    } else if (geoMultiscaleProfile == "uniform") {
      multiscaleProfile = AxialGeoMultiscaleMapping::MultiscaleProfile::UNIFORM;
    } else {
      PRECICE_UNREACHABLE("Unknown cross-section profile \"{}\".", geoMultiscaleProfile);
    }
 
    AxialGeoMultiscaleMapping::MultiscaleCrossSection multiscaleCrossSection = AxialGeoMultiscaleMapping::MultiscaleCrossSection::CIRCLE;
    if (geoMultiscaleCrossSection == "circle") {
      multiscaleCrossSection = AxialGeoMultiscaleMapping::MultiscaleCrossSection::CIRCLE;
    } else if (geoMultiscaleCrossSection == "square") {
      multiscaleCrossSection = AxialGeoMultiscaleMapping::MultiscaleCrossSection::SQUARE;
    } else {
      PRECICE_UNREACHABLE("Unknown geometric cross section \"{}\".", geoMultiscaleCrossSection);
    }
 
    configuredMapping.mapping = PtrMapping(new AxialGeoMultiscaleMapping(constraintValue, fromMesh->getDimensions(), multiscaleDimension, multiscaleType, multiscaleAxis, multiscaleRadius, multiscaleProfile, multiscaleCrossSection));
 
  } else if (type == TYPE_RADIAL_GEOMETRIC_MULTISCALE) {
 
    // Radial geometric multiscale is not applicable with the conservative constraint
    PRECICE_CHECK(constraintValue != Mapping::CONSERVATIVE,
                  "Radial geometric multiscale mapping is not implemented for the \"conservative\" constraint. "
                  "Please select constraint=\" consistent\" or a different mapping method.");
 
    // Convert strings into enums
    RadialGeoMultiscaleMapping::MultiscaleAxis multiscaleAxis;
    if (geoMultiscaleAxis == "x") {
      multiscaleAxis = RadialGeoMultiscaleMapping::MultiscaleAxis::X;
    } else if (geoMultiscaleAxis == "y") {
      multiscaleAxis = RadialGeoMultiscaleMapping::MultiscaleAxis::Y;
    } else if (geoMultiscaleAxis == "z") {
      multiscaleAxis = RadialGeoMultiscaleMapping::MultiscaleAxis::Z;
    } else {
      PRECICE_UNREACHABLE("Unknown geometric multiscale axis \"{}\".", geoMultiscaleAxis);
    }
 
    RadialGeoMultiscaleMapping::MultiscaleType multiscaleType;
    if (geoMultiscaleType == "spread") {
      multiscaleType = RadialGeoMultiscaleMapping::MultiscaleType::SPREAD;
    } else if (geoMultiscaleType == "collect") {
      multiscaleType = RadialGeoMultiscaleMapping::MultiscaleType::COLLECT;
    } else {
      PRECICE_UNREACHABLE("Unknown geometric multiscale type \"{}\".", geoMultiscaleType);
    }
 
    configuredMapping.mapping = PtrMapping(new RadialGeoMultiscaleMapping(constraintValue, fromMesh->getDimensions(), multiscaleType, multiscaleAxis));
 
  } else {
    // We need knowledge about the basis function in order to instantiate the rbf related mapping
    PRECICE_ASSERT(requiresBasisFunction(type));
    configuredMapping.mapping                = nullptr;
    configuredMapping.configuredWithAliasTag = type == TYPE_RBF_ALIAS;
  }
 
  configuredMapping.requiresBasisFunction = requiresBasisFunction(type);
 
  return configuredMapping;
}
 
void MappingConfiguration::checkDuplicates(const ConfiguredMapping &mapping)
{
  for (const ConfiguredMapping &configuredMapping : _mappings) {
    bool sameToMesh   = mapping.toMesh->getName() == configuredMapping.toMesh->getName();
    bool sameFromMesh = mapping.fromMesh->getName() == configuredMapping.fromMesh->getName();
    bool sameMapping  = sameToMesh && sameFromMesh;
    PRECICE_CHECK(!sameMapping,
                  "There cannot be two mappings from mesh \"{}\" to mesh \"{}\". "
                  "Please remove one of the duplicated meshes. ",
                  mapping.fromMesh->getName(), mapping.toMesh->getName());
  }
  for (const ConfiguredMapping &configuredMapping : _mappings) {
    bool sameToMesh  = mapping.toMesh->getName() == configuredMapping.toMesh->getName();
    bool isWrite     = mapping.direction == mapping::MappingConfiguration::WRITE && configuredMapping.direction == mapping::MappingConfiguration::WRITE;
    bool sameMapping = sameToMesh && isWrite && (mapping.fromMesh->isJustInTime() || configuredMapping.fromMesh->isJustInTime());
    PRECICE_CHECK(!sameMapping,
                  "There cannot be two mappings to mesh \"{}\". "
                  "Here, we have a mixture of just-in-time mapping and a conventional mapping. ",
                  mapping.toMesh->getName());
  }
}
 
void MappingConfiguration::xmlEndTagCallback(const xml::ConfigurationContext &context, xml::XMLTag &tag)
{
  if (tag.getNamespace() == TAG) {
    if (requiresBasisFunction(tag.getName())) {
      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());
      if (!_executorConfig) {
        _executorConfig = std::make_unique<ExecutorConfiguration>();
      }
      finishRBFConfiguration();
      _executorConfig.reset();
    }
    PRECICE_ASSERT(_mappings.back().mapping != nullptr);
  }
}
 
void MappingConfiguration::finishRBFConfiguration()
{
  PRECICE_ASSERT(_executorConfig);
  ConfiguredMapping &mapping = _mappings.back();
 
  // Instantiate the RBF mapping classes
  // We first categorize according to the executor
  // 1. the CPU executor
  if (_executorConfig->executor == ExecutorConfiguration::Executor::CPU) {
    if (_rbfConfig.solver == RBFConfiguration::SystemSolver::GlobalDirect) {
      mapping.mapping = getRBFMapping<RBFBackend::Eigen>(_rbfConfig.basisFunction, constraintValue, mapping.fromMesh->getDimensions(), _rbfConfig.supportRadius, _rbfConfig.shapeParameter, _rbfConfig.deadAxis, _rbfConfig.polynomial);
    } else if (_rbfConfig.solver == RBFConfiguration::SystemSolver::GlobalIterative) {
#ifndef PRECICE_NO_PETSC
      // for petsc initialization
      utils::Petsc::initialize(utils::Parallel::current()->comm);
 
      mapping.mapping = getRBFMapping<RBFBackend::PETSc>(_rbfConfig.basisFunction, constraintValue, mapping.fromMesh->getDimensions(), _rbfConfig.supportRadius, _rbfConfig.shapeParameter, _rbfConfig.deadAxis, _rbfConfig.solverRtol, _rbfConfig.polynomial);
#else
      PRECICE_ERROR("The global-iterative RBF solver on a CPU requires a preCICE build with PETSc enabled.");
#endif
    } else if (_rbfConfig.solver == RBFConfiguration::SystemSolver::PUMDirect) {
      _ginkgoParameter          = GinkgoParameter();
      _ginkgoParameter.executor = "cpu";
      mapping.mapping           = getRBFMapping<RBFBackend::PUM>(_rbfConfig.basisFunction, constraintValue, mapping.fromMesh->getDimensions(), _rbfConfig.supportRadius, _rbfConfig.shapeParameter, _rbfConfig.polynomial, _rbfConfig.verticesPerCluster, _rbfConfig.relativeOverlap, _rbfConfig.projectToInput, _ginkgoParameter);
    } else {
      PRECICE_UNREACHABLE("Unknown RBF solver.");
    }
    // 2. any other executor is configured via Ginkgo
  } else {
    _ginkgoParameter                   = GinkgoParameter();
    _ginkgoParameter.usePreconditioner = false;
    _ginkgoParameter.deviceId          = _executorConfig->deviceId;
    if (_executorConfig->executor == ExecutorConfiguration::Executor::CUDA) {
      _ginkgoParameter.executor = "cuda-executor";
#ifndef PRECICE_WITH_CUDA
      PRECICE_ERROR("The cuda-executor (configured for the mapping from mesh {} to mesh {}) requires a Kokkos and preCICE build with Cuda enabled.", mapping.fromMesh->getName(), mapping.toMesh->getName());
#endif
    } else if (_executorConfig->executor == ExecutorConfiguration::Executor::HIP) {
      _ginkgoParameter.executor = "hip-executor";
#ifndef PRECICE_WITH_HIP
      PRECICE_ERROR("The hip-executor (configured for the mapping from mesh {} to mesh {}) requires a Kokkos and preCICE build with HIP enabled.", mapping.fromMesh->getName(), mapping.toMesh->getName());
#endif
    } else if (_executorConfig->executor == ExecutorConfiguration::Executor::SYCL) {
      _ginkgoParameter.executor = "sycl-executor";
#ifndef PRECICE_WITH_SYCL
      PRECICE_ERROR("The sycl-executor (configured for the mapping from mesh {} to mesh {}) requires a Kokkos and preCICE build with SYCL enabled.", mapping.fromMesh->getName(), mapping.toMesh->getName());
#endif
    } else if (_executorConfig->executor == ExecutorConfiguration::Executor::OpenMP) {
      _ginkgoParameter.executor = "omp-executor";
      _ginkgoParameter.nThreads = _executorConfig->nThreads;
#ifndef PRECICE_WITH_OPENMP
      PRECICE_ERROR("The omp-executor (configured for the mapping from mesh {} to mesh {}) requires a Kokkos and preCICE build with OpenMP enabled.", mapping.fromMesh->getName(), mapping.toMesh->getName());
#endif
    }
    if (_rbfConfig.solver == RBFConfiguration::SystemSolver::GlobalDirect) {
#ifndef PRECICE_NO_GINKGO
      _ginkgoParameter.solver = "qr-solver";
      mapping.mapping         = getRBFMapping<RBFBackend::Ginkgo>(_rbfConfig.basisFunction, constraintValue, mapping.fromMesh->getDimensions(), _rbfConfig.supportRadius, _rbfConfig.shapeParameter, _rbfConfig.deadAxis, _rbfConfig.polynomial, _ginkgoParameter);
#else
      PRECICE_ERROR("The selected direct solver for the global RBF mapping on executor {} from mesh {} to mesh {} requires a preCICE build with Ginkgo enabled.", _ginkgoParameter.executor, mapping.fromMesh->getName(), mapping.toMesh->getName());
#endif
    } else if (_rbfConfig.solver == RBFConfiguration::SystemSolver::GlobalIterative) {
#ifndef PRECICE_NO_GINKGO
      _ginkgoParameter.solver       = "cg-solver";
      _ginkgoParameter.residualNorm = _rbfConfig.solverRtol;
      mapping.mapping               = getRBFMapping<RBFBackend::Ginkgo>(_rbfConfig.basisFunction, constraintValue, mapping.fromMesh->getDimensions(), _rbfConfig.supportRadius, _rbfConfig.shapeParameter, _rbfConfig.deadAxis, _rbfConfig.polynomial, _ginkgoParameter);
#else
      PRECICE_ERROR("The selected iterative solver for the global RBF mapping on executor {} from mesh {} to mesh {} requires a preCICE build with Ginkgo enabled.", _ginkgoParameter.executor, mapping.fromMesh->getName(), mapping.toMesh->getName());
#endif
    } else if (_rbfConfig.solver == RBFConfiguration::SystemSolver::PUMDirect) {
#ifndef PRECICE_NO_KOKKOS_KERNELS
      PRECICE_CHECK(!(mapping.fromMesh->isJustInTime() || mapping.toMesh->isJustInTime()), "Executor \"{}\" is not implemented as just-in-time mapping.", _ginkgoParameter.executor);
      mapping.mapping = getRBFMapping<RBFBackend::PUM>(_rbfConfig.basisFunction, constraintValue, mapping.fromMesh->getDimensions(), _rbfConfig.supportRadius, _rbfConfig.shapeParameter, _rbfConfig.polynomial, _rbfConfig.verticesPerCluster, _rbfConfig.relativeOverlap, _rbfConfig.projectToInput, _ginkgoParameter, _executorConfig->computeEvaluationOffline);
#else
      PRECICE_ERROR("The selected pu-rbf solver using executor \"{}\" for the mapping from mesh {} to mesh {} requires a preCICE build with Kokkos-kernels enabled.", _ginkgoParameter.executor, mapping.fromMesh->getName(), mapping.toMesh->getName());
#endif
    } else {
      PRECICE_UNREACHABLE("Unknown solver type.");
    }
  }
}
 
const std::vector<MappingConfiguration::ConfiguredMapping> &MappingConfiguration::mappings()
{
  return _mappings;
}
 
bool MappingConfiguration::requiresBasisFunction(const std::string &mappingType) const
{
  return mappingType == TYPE_RBF_PUM_DIRECT || mappingType == TYPE_RBF_GLOBAL_DIRECT || mappingType == TYPE_RBF_GLOBAL_ITERATIVE || mappingType == TYPE_RBF_ALIAS;
}
 
BasisFunction MappingConfiguration::parseBasisFunctions(const std::string &basisFctName) const
{
  BasisFunction basisFunction{};
  if (basisFctName == RBF_TPS)
    basisFunction = BasisFunction::ThinPlateSplines;
  else if (basisFctName == RBF_MULTIQUADRICS)
    basisFunction = BasisFunction::Multiquadrics;
  else if (basisFctName == RBF_INV_MULTIQUADRICS)
    basisFunction = BasisFunction::InverseMultiquadrics;
  else if (basisFctName == RBF_VOLUME_SPLINES)
    basisFunction = BasisFunction::VolumeSplines;
  else if (basisFctName == RBF_GAUSSIAN)
    basisFunction = BasisFunction::Gaussian;
  else if (basisFctName == RBF_CTPS_C2)
    basisFunction = BasisFunction::CompactThinPlateSplinesC2;
  else if (basisFctName == RBF_CPOLYNOMIAL_C0)
    basisFunction = BasisFunction::WendlandC0;
  else if (basisFctName == RBF_CPOLYNOMIAL_C2)
    basisFunction = BasisFunction::WendlandC2;
  else if (basisFctName == RBF_CPOLYNOMIAL_C4)
    basisFunction = BasisFunction::WendlandC4;
  else if (basisFctName == RBF_CPOLYNOMIAL_C6)
    basisFunction = BasisFunction::WendlandC6;
  else if (basisFctName == RBF_CPOLYNOMIAL_C8)
    basisFunction = BasisFunction::WendlandC8;
  else
    PRECICE_UNREACHABLE("Unknown basis function \"{}\".", basisFctName);
  return basisFunction;
}
} // namespace precice::mapping