Build your own adapter

The deal.II adapter provides a minimal example of a deal.II code, which has been coupled using preCICE. Since most users will probably have their own deal.II-based solvers, this section explains the preCICE-related code changes. A more general adapter example with a step-by-step tutorial is also available on the preCICE wiki.

Contact us if you have any questions. Even if you don't have any questions, please let us know about your experience when your adapter is ready!

Which information is needed by preCICE?

preCICE uses a black-box coupling approach, which means the solver only needs to provide a minimal set of information. In the simplest case, this includes configuration information e.g. name of the participant and the coordinates of the data you want to exchange, e.g. the mesh vertices. If you want to use a nearest-projection mapping, you need to additionally specify mesh connectivity between the vertices, which is currently not included in this adapter example.

About fluid-structure coupling

For every multi-physics coupling, proper coupling data needs to be exchanged between all participants. In our example case, the Fluid participant calculates forces (per cell face), which are passed to the Solid participant. Using the forces for the structural calculations, the Solid participant calculates displacements, which are then passed back to the Fluid participant. As outlined above, preCICE needs coordinates of the data points you want to exchange. Since the forces live on the cell faces and the displacements live on the cell vertices, we define two coupling meshes in this example: A face-based mesh for the forces: Fluid -> Solid and a vertex based mesh for the displacements: Solid -> Fluid. This is not necessarily the case and depends on your spatial discretization method and the coupling data. If both data sets live on the same mesh element, defining one mesh would be sufficient.

Step-by-step guidance

This guideline can be used to couple your own code using preCICE. All code snippets have been taken from the coupled_elasto_dynamics.cc file.

Step 1: Prepare your solver

As a first preparation step, you need to include the preCICE library headers, which is simply the SolverInterface.hpp

#include "precice/SolverInterface.hpp"

The preCICE API is explained in the source code.

Step 2: Set up API

The API is located in precice::SolverInterface, and its constructor requires the participant's name, and the rank and size of the current thread. Since it is typical for deal.II to create a class for each code and call the run() function so start the simulation, the precice::SolverInterface object is generated in the constructor of the CoupledElastoDynamics class. Once the basic preCICE interface is set up, we can steer the behaviour of preCICE.

Step 3: Initialize preCICE

Initializing preCICE is done in the respective initialize_precice() function, which is located in the run() function. Here, the configuration file is read and the number of interface nodes and the interface coordinates are passed to preCICE. Have a look in the source code and the comments for the implementation. Depending on the location of your coupling data (cf. fluid-structure section), you might be able to copy this function completely.

Note 1: In order to exchange data between participants, preCICE needs to know which data is actually in the data array. Therefore, the data names in the precice-config.xml are converted to preCICE specific IDs. Make sure the naming in the respective functions is the same in the source code and in the configuration file. For example, if you name the coupling data force Force in your precice-config.xml, you need to pass the same name Force in the getDataID(Force,...) function.

Note 2: Several declarations for preCICE related functions or variables have been made in the constructor of the main class.

Note 3: Apart from the coordinates of the data points, we need to pass the number of data points to preCICE. In case of the interface faces, this step is included at the end of the make_grid() function, where the boundary_IDs are set as well, in order to avoid code duplication.

Note 4: For each information you want to pass to preCICE, you need to provide the information in a specific data format. Therefore, all vector components need to be listed in an array as shown in the following example. Think of the data coordinates: Each coordinate is a dim dimensional vector with v1 = [x1; y1; z1] , v2 = [x2; y2; z2] ... For preCICE, all components are collected in a single array and arranged like: coordinates = [x1 y1 z1 x2 y2 z2 x3...].

Note 5: We need to use the same order for our coupling data, in which we pass the coordinates to preCICE. Further information are mentioned in the next step.

Note 6: On several code sections, you will find an if-statement for the enable precice option. This has been introduced to allow the user preliminary computations without preCICE. A constant value can be set to compute simple testcases and e.g. have a look at the grid. Feel free to simply ignore this for your own adapter.

Step 4: Modify the time loop

Coupled problems are usually time dependent. We set up a separate Time class, which handles all time related tasks. preCICE uses specific functions in the time loop, which steer the simulation. The time loop is located in the compute_timesteps() function. The relevant loop for our purpose is the enable_precice = true loop. `Therefore, have a look at the following simplified basic version of our example program:

     while(precice.isCouplingOngoing() &&
           time.get_timestep() < time.get_n_timesteps())  // ask preCICE for the coupling process
     {
         time.increment();       // increment time 

         assemble_rhs();         // apply coupling data in the rhs 

         solve();                // solve the system as ususal

         update_displacement();  // update time dependent variables as usual

         advance_precice();      // exchange coupling data

         output_results(time.get_timestep()); // store result files as usual
     }

This time loop would be suited for an explicit coupling.

assemble_rhs(): Since now the coupling data is time dependent, you need to assemble parts of your matrices/vectors in every time step, depending on your coupling problem/implementation. In most of the cases and in our example, this will be the rhs of the equations. We obtain coupling data (depending on your preCICE configuration) before the first step and at the end of each time step. The function assemble_rhs() uses the coupling data, namely forces per cell face, to rebuild the rhs vector. How to apply coupling data is strongly problem-dependent. In our case, we use a face quadrature to apply a Neumann boundary condition. Have a look in the source code for the implementation. In addition to Note 5 (previous step), make sure to stay consistent between the extraction of the face coordinates and the reading of the force data (named precice_forces in the source code). In the example this is ensured by iterating in the same way over all cells by using the same cell_iterator().

advance_precice(): After we computed a time step, we want to exchange data with other participants. This is done in the advance_precice() function. In the beginning we provide our coupling data, displacements in the example, and obtain coupling data after calling precice.advance(time.get_delta_t()). This function will be the same in most of the adapters. But in order to provide your coupling data, you need to extract the coupling data from your global solution vector - again in a consistent order with the data coordinates you passed to preCICE in the initialize_precice() function. We wrote an additional function called extract_relevant_displacements for this task. If you have the same data structure for your global solution vector, you should be able to use a similar approach for your project.

Step 5: Terminating preCICE

After all time steps have been computed, the function finalize() is used to free relevant data structures and close the communication channel.

Step 6: Support implicit coupling

The steps 1-5 are sufficient for an explicit coupling. In many cases, an implicit coupling is desired. Therefore, we need to save the time-dependent data before we start recomputing it and reload it later, if we have not yet reached a certain convergence criterion. This leads to two additional functions wrapped around an if-statement in the time loop:

     while(precice.isCouplingOngoing() &&
           time.get_timestep() < time.get_n_timesteps())
     {
         time.increment();

         if (precice.isActionRequired(precice::constants::actionWriteIterationCheckpoint())){
             save_old_state();
             precice.markActionFulfilled(precice::constants::actionWriteIterationCheckpoint());
         }

         assemble_rhs();

         solve();

         update_displacement();

         advance_precice();

         if (precice.isActionRequired(precice::constants::actionReadIterationCheckpoint())){
             reload_old_state();
             precice.markActionFulfilled(precice::constants::actionReadIterationCheckpoint());
         }

         if(precice.isTimeWindowComplete())
             output_results(time.get_timestep());
     }

The two functions save_old_state() and reload_old_state() copy the time dependent variables and reload them later. Note, that we don't simply use the 'old' variables to reload the older state, since we additionally want to support subcycling, where one participant e.g. solid, computes multiple solver time steps during a coupling time step. The time loop structure for implicit couplings will probably be the same in your own project.

Your adapter is now ready for most of the preCICE features. For nearest-projection mapping, mesh connectivity needs to be provided as well. An idea of how this could be implemented will be soon available in our issues.