[MappingApplication] Creation of a Radial Basis Function Mapper for IBRA and FEM

[juancamarotti]

This PR extends the MappingApplication to support Radial Basis Function (RBF) mapping, enabling flexible, mesh-independent transfer of scalar and vector fields between non-matching geometries.

Overview

• Introduces an RBFMapper for scalar and vector field transfer across arbitrary point sets, surfaces, and volumes.
• Supports multiple RBF kernels (Gaussian, multiquadric, inverse multiquadric, thin-plate spline) with optional polynomial augmentation.
• Implements bi-directional mapping using consistent interpolation matrices. (Note: Consistent mapping is supported from IBRA → FEM domains or between FEM domains.)
• Fully integrated with the CoSimulationApplication for FSI and other multi-physics couplings.
• Supports conservative mapping, ensuring energy and load consistency.

Tests

• Validated with analytical scalar and vector mapping benchmarks.
• Verified on FSI benchmarks, including FSI Mok and FSI Turek, demonstrating accuracy and robustness.

Multivariate interpolation for fluid-structure-interaction problems using RBF.pdf

[philbucher]

for now I did only a brief review.
What confuses me most is that you are never solving a linear system, which to my understanding is required for the rbfs, how else would you compute the coefficients?
Thats one of the main reasons I never implemented it myself

[philbucher]

mapping is physics-agnosic, pls change the naming

[juancamarotti]

Changed 👍

[juancamarotti]

Hi @philbucher, I will try to clarify your point regarding the solution of the linear system involved by the RBF mapper.

This mapper uses Radial Basis Functions (RBFs) to build a linear mapping between origin and destination nodes.
The idea is simple:

1. Build the RBF interpolation matrix between the origin points --> Φorigin​[i,j]=ϕ(∥xiorigin​−xjorigin​∥)
2. Invert that matrix once at the beginning of the simulation
3. Build another RBF matrix that relates destination points to the origin points --> Φdest​[i,j]=ϕ(∥xidest​−xjorigin​∥)
4. Multiply them to get a precomputed mapping matrix --> M=Φdest​ * [Φorigin​]−1
5. Then, at every time step or mapping call, you just do --> y_dest​=M * y_origin​
   So the heavy part (the inversion) happens only once during setup, and applying the mapping later is just a fast matrix–vector multiplication.

So the heavy part (the inversion) happens only once during setup, and applying the mapping later is just a fast matrix–vector multiplication.

Of course, there are other mapping techniques that can be less computationally expensive. However, it’s still valuable to have this RBF mapper available in Kratos as it’s a standard approach widely used in low-fidelity aeroelastic coupling frameworks in the industry (for example, in Airbus).

The RBF system can’t be solved ahead of time because it would require knowing the origin displacements, which are only available during the simulation. Instead, we precompute a linear mapping matrix using just the origin and destination geometries. At each time step, the mapping is applied with a simple matrix–vector multiplication, avoiding the need to solve a system repeatedly.

It's worth commenting that usually in low-fidelity aeroelastic frameworks, the size of this matrix to be inverted is not so big, bacause the loads are not applied in every node of the structure but just in some specific stiff nodes.

[philbucher]

thanks for the explanation

So the heavy part (the inversion) happens only once during setup, and applying the mapping later is just a fast matrix–vector multiplication.

Oh boy 😅

for tiny systems it might work, but it will blow up very soon. I will write a longer reply with my thoughts later

[juancamarotti]

thanks for the explanation

So the heavy part (the inversion) happens only once during setup, and applying the mapping later is just a fast matrix–vector multiplication.

Oh boy 😅

for tiny systems it might work, but it will blow up very soon. I will write a longer reply with my thoughts later

I know this is non optimal, but I really doubt someone would use this mapper for a large-scale simulation... if you have a better solution for this, I will be very happy to hear, and I still think this mapper is a good addition to the MappingApplication. Thanks for the help!

[philbucher]

After reviewing this and giving it some thought I see two ways forward:

• You pretty much merge as is, but rename it to sth like BacisRbfMapper or SimpleRbfMapper. I am quite convinced that precomputing the mapping-matrix via inversion is not a good approach (and it doesnt work in MPI). Since I am not in favor of this and no longer a user of Kratos, I would ask @sunethwarna or someone from the tech-committee for approval.
  Whats the largest example you tried, what was the size of the mapping-matrix? Have you done benchmarking or profiling? I suggest you try with sizes of examples you want to use in your thesis, to not realize that it doesnt work while running the last case for results ;)

• Alternatively, I could help you in implementing this properly and adhering to the existing design of the MappingApp. This would be my recommended approach, as it would also support MPI and be much better in terms of performance, also for big problems. It would look very similar to the barycentric mapper, plus additionally a linear solution in every mapping step.
  Even with this approach you can precompute the matrix, check the CouplingGeometriesMapper, this one supports both options. There we first tried pretty much your approach but realized that it was not working for real problems.

Up to you, let me know :)

BTW the RBF mapper is default in precice, where they use it for huuuuge examples.

[juancamarotti]

After reviewing this and giving it some thought I see two ways forward:

• You pretty much merge as is, but rename it to sth like BacisRbfMapper or SimpleRbfMapper. I am quite convinced that precomputing the mapping-matrix via inversion is not a good approach (and it doesnt work in MPI). Since I am not in favor of this and no longer a user of Kratos, I would ask @sunethwarna or someone from the tech-committee for approval.
  Whats the largest example you tried, what was the size of the mapping-matrix? Have you done benchmarking or profiling? I suggest you try with sizes of examples you want to use in your thesis, to not realize that it doesnt work while running the last case for results ;)
• Alternatively, I could help you in implementing this properly and adhering to the existing design of the MappingApp. This would be my recommended approach, as it would also support MPI and be much better in terms of performance, also for big problems. It would look very similar to the barycentric mapper, plus additionally a linear solution in every mapping step.
  Even with this approach you can precompute the matrix, check the CouplingGeometriesMapper, this one supports both options. There we first tried pretty much your approach but realized that it was not working for real problems.

Up to you, let me know :)

BTW the RBF mapper is default in precice, where they use it for huuuuge examples.

Thank you for the detailed feedback @philbucher! I agree that the second approach is the cleaner and more robust solution, especially considering MPI support and performance for large problems.

So far, the only benchmarking I’ve done was with the FSI Mok and Turek examples. They’re good for testing but relatively small, so I haven’t fully assessed performance for larger cases yet.

It would be great to have your guidance and help for implementing this properly. What do you think?

[philbucher]

sounds good!

So there are a few things:

• How big do want the search-radius to be? I.e. how global should the support of the RBFs be? => this decision will impact how the search-settings will be specified, i.e. slightly similar to the one of the other sovlers
• How many points do you want to consider for the RBFs? => the more you consider, the denser the matrix will be
• You basically collect the points you want via the InterfaceInfo, and then compute the RBFs with a LocalSystem. Conceptually the Barycentric mapper is very similar, check for example BarycentricLocalSystem::CalculateAll
• in more detail:
  • Define rbf-support-radius (aka search-radius) and num support points, and pass it to the RbfInterfaceInfo. See BarycentricInterfaceInfo::ProcessSearchResult. The search structures handle the rest, both in serial and MPI
  • in RbfLocalSystem you first combine the results from all RbfInterfaceInfo (one from each rank, basically the same as is done in ProcessSearchResult:

    Kratos/applications/MappingApplication/custom_mappers/barycentric_mapper.cpp

    Lines 276 to 280 in 601a4fb

    	for (std::size_t i=1; i<mInterfaceInfos.size(); ++i) {
    	// in MPI there might be results also from other ranks
    	const BarycentricInterfaceInfo& r_info = static_cast<const BarycentricInterfaceInfo&>(*mInterfaceInfos[i]);
    	closest_points.Merge(r_info.GetClosestPoints());
    	}

  • Then you have the points you need to compute the RBFs. Insert them in the local mapping matrix, then they will be assembled into the coefficients matrix by the MappingMatrixBuilder (again, it takes care of the MPI stuff)
  • Last step is to either precompute the matrix, or add a linear every time you map. Same as is done for the CouplingGeometryMapper:

    Kratos/applications/MappingApplication/custom_mappers/coupling_geometry_mapper.cpp

    Lines 230 to 256 in 601a4fb

    	template<class TSparseSpace, class TDenseSpace>
    	void CouplingGeometryMapper<TSparseSpace, TDenseSpace>::MapInternal(
    	const Variable<double>& rOriginVariable,
    	const Variable<double>& rDestinationVariable,
    	Kratos::Flags MappingOptions)
    	{
    	const bool dual_mortar = mMapperSettings["dual_mortar"].GetBool();
    	const bool precompute_mapping_matrix = mMapperSettings["precompute_mapping_matrix"].GetBool();
    	mpInterfaceVectorContainerMaster->UpdateSystemVectorFromModelPart(rOriginVariable, MappingOptions);
    	if (dual_mortar || precompute_mapping_matrix) {
    	TSparseSpace::Mult(
    	*mpMappingMatrix,
    	mpInterfaceVectorContainerMaster->GetVector(),
    	mpInterfaceVectorContainerSlave->GetVector()); // rQd = rMdo * rQo
    	} else {
    	TSparseSpace::Mult(
    	*mpMappingMatrixProjector,
    	mpInterfaceVectorContainerMaster->GetVector(),
    	*mpTempVector); // rQd = rMdo * rQo
    	mpLinearSolver->Solve(*mpMappingMatrixSlave, mpInterfaceVectorContainerSlave->GetVector(), *mpTempVector);
    	}
    	mpInterfaceVectorContainerSlave->UpdateModelPartFromSystemVector(rDestinationVariable, MappingOptions);
    	}

That should be it ...

In the meantime can you make a separate PR for the changes in the RBF utility?

[philbucher]

this MR is outdated and is not planned to be merged, right?

[juancamarotti]

this MR is outdated and is not planned to be merged, right?

Hi @philbucher. Yes, this PR is only a draft and thus is not planned to be merged. The one to be merged is the following PR #13981.