deal.II version GIT relicensing-2287-g6548a49e0a 2024-12-20 18:30:00+00:00
|
This program was contributed by David Schneider <dav.schneider@tum.de>.
It comes without any warranty or support by its authors or the authors of deal.II.
This program is part of the deal.II code gallery and consists of the following files (click to inspect):
preCICE allows to couple deal.II to external simulation software, such as OpenFOAM, SU2, or CalculiX. To keep dependencies of this example minimal, we couple deal.II to an external C++ program, which provides a time varying boundary condition. The deal.II code consists mainly of the step-4 tutorial program, where a simple Laplace problem is solved.
Coupling with preCICE is usually carried out along surfaces in order to apply a Dirichlet-Neumann coupling between two domains (volume coupling is also possible). For the sake of simplicity, we couple here an external C++ program in a unidirectional fashion to one side of our quadrilateral domain. The external C++ program generates a parabolic boundary profile with time varying amplitude. The boundary values are then used in the Laplace solver as a Dirichlet boundary condition.
Coupled simulations deal mostly with time-dependent problems. Hence, we make the stationary Laplace problem from step-4 time dependent.
\begin{align*} \frac{\partial u}{\partial t}-\Delta u &= f \qquad\qquad & \text{in}\ \Omega, \\ u &= x^2+y^2 \qquad\qquad & \text{on}\ \partial\Omega_s, \\ u &= g(t) \qquad\qquad & \text{on}\ \partial\Omega_c. \end{align*}
with the fixed Dirichlet boundary \Omega_s, the coupling boundary \Omega_c and the time-dependent coupling data g(t).
The system is consequently discretized by a first-order backward Euler method, resulting in
\begin{align*} \frac{u^{n+1}-u^n}{\Delta t}-\Delta u^{n+1} &= f \qquad\qquad & \text{in}\ \Omega, \end{align*}
at the next time level n+1.
deal.II
, version 9.2
or greater. Older versions might work as well, but have not been tested.
preCICE, version 2.0
or greater. Have a look at the provided link for an installation guide.
Similar to the example programs, run
in this directory to configure the problem. The explicit specification of the library locations can be omitted if they are installed globally or via environment variables. You can then switch between debug and release mode by calling either
or
This command will generate two executables: one for the coupled_laplace_problem
and one for the fancy_boundary_condition
participant.
executes the coupled_laplace_problem
. In order to start the coupled simulation, execute
in the same directory from another terminal window.
preCICE and the deal.II solver create several log files during the simulation. In order to remove all result-associated files and clean-up the simulation directory the Allclean script can be executed.
By default, the results are saved in each time step using vtk
files, which can, for example, be visualized with ParaView. You should get something like:
A complete overview of the preCICE project can be found on the preCICE webpage. The source code of preCICE is hosted on Github. For more coupled deal.II codes, have a look at the deal.II adapter repository. In the preCICE tutorials, a fluid-structure interaction example is given coupling OpenFOAM (fluid dynamics) to deal.II (solid mechanics). The preCICE reference paper
The included deal.II header files are the same as in the other example programs:
In addition to the deal.II header files, we include the preCICE API in order to obtain access to preCICE specific functionality
Configuration parameters
We set up a simple hard-coded struct containing all names we need for external coupling. The struct includes the name of the preCICE configuration file as well as the name of the simulation participant, the name of the coupling mesh and the name of the exchanged data. The last three names you also find in the preCICE configuration file. For real application cases, these names are better handled by a parameter file.
The Adapter class
The Adapter class handles all functionalities to couple the deal.II solver code to other solvers with preCICE, i.e., data structures are set up and all relevant information is passed to preCICE.
public precCICE solver interface
Boundary ID of the deal.II triangulation, associated with the coupling interface. The variable is defined in the constructor of this class and intentionally public so that it can be used during the grid generation and system assembly. The only thing, one needs to make sure is that this ID is unique for a particular triangulation.
preCICE related initializations These variables are specified in and read from a parameter file, which is in this simple tutorial program the CouplingParameter struct already introduced in the beginning.
The node IDs are filled by preCICE during the initialization and associated to the spherical vertices we pass to preCICE
DoF IndexSet, containing relevant coupling DoF indices at the coupling boundary
Data containers which are passed to preCICE in an appropriate preCICE specific format
The MPI rank and total number of MPI ranks is required by preCICE when the Participant is created. Since this tutorial runs only in serial mode we define the variables manually in this class instead of using the regular MPI interface.
Function to transform the obtained data from preCICE into an appropriate map for Dirichlet boundary conditions
In the constructor of the Adapter class, we set up the preCICE Participant. We need to tell preCICE our name as participant of the simulation and the name of the preCICE configuration file. Both have already been specified in the CouplingParameter class above. Thus, we pass the class directly to the constructor and read out all relevant information. As a second parameter, we need to specify the boundary ID of our triangulation, which is associated with the coupling interface.
This function initializes preCICE (e.g. establishes communication channels and allocates memory) and passes all relevant data to preCICE. For surface coupling, relevant data is in particular the location of the data points at the associated interface(s). The boundary_data
is an empty map, which is filled by preCICE, i.e., information of the other participant. Throughout the system assembly, the map can directly be used in order to apply the Dirichlet boundary conditions in the linear system.
Afterwards, we extract the number of interface nodes and the coupling DoFs at the coupling interface from our deal.II solver via extract_boundary_dofs()
The ComponentMask()
might be important in case we deal with vector valued problems, because vector valued problems have a DoF for each component.
The coupling DoFs are used to set up the boundary_data
map. At the end, we associate here each DoF with a respective boundary value.
Since we deal with a scalar problem, the number of DoFs at the particular interface corresponds to the number of interface nodes.
Now, we need to tell preCICE the coordinates of the interface nodes. Hence, we set up a std::vector to pass the node positions to preCICE. Each node is specified only once.
Set up the appropriate size of the data container needed for data exchange. Here, we deal with a scalar problem, so that only a scalar value is read/written per interface node.
The IDs are filled by preCICE during the initializations.
The node location is obtained using map_dofs_to_support_points()
.
support_points
contains now the coordinates of all DoFs. In the next step, the relevant coordinates are extracted using the IndexSet with the extracted coupling_dofs.
Now we have all information to define the coupling mesh and pass the information to preCICE.
Then, we initialize preCICE internally calling the API function initialize()
here, we obtain data, i.e. the boundary condition, from another participant. We have already vertex IDs and just need to convert our obtained data to the deal.II compatible 'boundary map' , which is done in the format_deal_to_precice function.
After receiving the coupling data in read_data_buffer
, we convert it to the std::map boundary_data
which is later needed in order to apply Dirichlet boundary conditions
The function advance()
is called in the main time loop after the computation in each time step. Here, preCICE exchanges the coupling data internally and computes mappings as well as acceleration methods.
We specify the computed time-step length and pass it to preCICE.
This function takes the std::vector obtained by preCICE in read_data_buffer
and inserts the values to the right position in the boundary map used throughout our deal.II solver for Dirichlet boundary conditions. The function is only used internally in the Adapter class and not called in the solver itself. The order, in which preCICE sorts the data in the read_data_buffer
vector is exactly the same as the order of the initially passed vertices coordinates.
We already stored the coupling DoF indices in the boundary_data
map, so that we can simply iterate over all keys in the map.
The solver class is essentially the same as in step-4. We only extend the stationary problem to a time-dependent problem and introduced the coupling. Comments are added at any point, where the workflow differs from step-4.
We allocate all structures required for the preCICE coupling: The map is used to apply Dirichlet boundary conditions and filled in the Adapter class with data from the other participant. The CouplingParameters hold the preCICE configuration as described above. The interface boundary ID is the ID associated to our coupling interface and needs to be specified, when we set up the Adapter class object, because we pass it directly to the Constructor of this class.
The time-step size delta_t is the actual time-step size used for all computations. The preCICE time-step size is obtained by preCICE in order to ensure a synchronization at all coupling time steps. The solver time step-size is the desired time-step size of our individual solver. In more sophisticated computations, it might be determined adaptively. The time_step
counter is just used for the time-step number.
We choose the boundary in positive x direction for the interface coupling.
Reset global structures
Update old solution values
The solution values from previous time steps are stored for each quadrature point
Get the local values from the ‘fe_values’ object
The system matrix contains additionally a mass matrix due to the time discretization. The RHS has contributions from the old solution values.
Copy local to global
At first, we apply the Dirichlet boundary condition from step-4, as usual.
Afterwards, we apply the coupling boundary condition. The boundary_data
has already been filled by preCICE.
After we set up the system, we initialize preCICE and the adapter using the functionalities of the Adapter.
preCICE steers the coupled simulation: isCouplingOngoing
is used to synchronize the end of the simulation with the coupling partner
The time step number is solely used to generate unique output files
preCICE returns the maximum admissible time-step size, which needs to be compared to our desired solver time-step size.
Next we read data. Since we use a fully backward Euler method, we want the data to be associated to the end of the current time-step (delta_t) Time-interpolation methods in preCICE allow to get readData at any point in time, if the coupling scheme allows it
In the time loop, we assemble the coupled system and solve it as usual.
After we solved the system, we advance the coupling to the next time level. In a bi-directional coupled simulation, we would pass our calculated data to preCICE
Write an output file if the time step is completed. In case of an implicit coupling, where individual time steps are computed more than once, the function isTimeWindowCompleted
prevents unnecessary result writing. For this simple tutorial configuration (explicit coupling), the function returns always true
.
This program does not use any deal.II functionality and depends only on preCICE and the standard libraries.
The program computes a time-varying parabolic boundary condition, which is passed to preCICE and serves as Dirichlet boundary condition for the other coupling participant.
Function to generate boundary values in each time step
Scale the current time value
Define the amplitude. Values run from -0.5 to 0.5
Specify the actual data we want to pass to the other participant. Here, we choose a parabola with boundary values 2 in order to enforce continuity to adjacent boundaries.
Configuration
Adjust to MPI rank and size for parallel computation
Set up data structures
Define a boundary mesh
The x-coordinate is always 1, i.e., the boundary is parallel to the y-axis. The y-coordinate is descending from 1 to -1.
Pass the vertices to preCICE
Variables for the time
Not used in the configuration by default
initialize the Participant
Start time loop
Generate new boundary data