--- /dev/null
+// The included deal.II header files are the same as in the other example
+// programs:
+#include <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+// In addition to the deal.II header files, we include the preCICE API in order
+// to obtain access to preCICE specific functionality
+#include <precice/SolverInterface.hpp>
+
+#include <fstream>
+#include <iostream>
+
+using namespace dealii;
+
+// 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.
+struct CouplingParamters
+{
+ const std::string config_file = "precice-config.xml";
+ const std::string participant_name = "laplace-solver";
+ const std::string mesh_name = "dealii-mesh";
+ const std::string read_data_name = "boundary-data";
+};
+
+
+// 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.
+
+template <int dim, typename ParameterClass>
+class Adapter
+{
+public:
+ Adapter(const ParameterClass & parameters,
+ const types::boundary_id dealii_boundary_interface_id);
+
+ double
+ initialize(const DoFHandler<dim> & dof_handler,
+ std::map<types::global_dof_index, double> &boundary_data,
+ const MappingQGeneric<dim> & mapping);
+
+ double
+ advance(std::map<types::global_dof_index, double> &boundary_data,
+ const double computed_timestep_length);
+
+ // public precCICE solver interface
+ precice::SolverInterface precice;
+
+ // 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.
+ const unsigned int dealii_boundary_interface_id;
+
+private:
+ // 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.
+ const std::string mesh_name;
+ const std::string read_data_name;
+
+ // These IDs are filled by preCICE during the initialization. We set a default
+ // value of -1 in order to detect potential errors more easily.
+ int mesh_id;
+ int read_data_id;
+ int n_interface_nodes;
+
+ // DoF IndexSet, containing relevant coupling DoF indices at the coupling
+ // boundary
+ IndexSet coupling_dofs;
+
+ // Data containers which are passed to preCICE in an appropriate preCICE
+ // specific format
+ std::vector<int> interface_nodes_ids;
+ std::vector<double> read_data;
+
+ // The MPI rank and total number of MPI ranks is required by preCICE when the
+ // SolverInterface 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.
+ static constexpr int this_mpi_process = 0;
+ static constexpr int n_mpi_processes = 1;
+
+ // Function to transform the obtained data from preCICE into an appropriate
+ // map for Dirichlet boundary conditions
+ void
+ format_precice_to_dealii(
+ std::map<types::global_dof_index, double> &boundary_data) const;
+};
+
+
+
+// In the constructor of the Adapter class, we set up the preCICE
+// SolverInterface. 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.
+template <int dim, typename ParameterClass>
+Adapter<dim, ParameterClass>::Adapter(
+ const ParameterClass & parameters,
+ const types::boundary_id deal_boundary_interface_id)
+ : precice(parameters.participant_name,
+ parameters.config_file,
+ this_mpi_process,
+ n_mpi_processes)
+ , dealii_boundary_interface_id(deal_boundary_interface_id)
+ , mesh_name(parameters.mesh_name)
+ , read_data_name(parameters.read_data_name)
+{}
+
+
+
+// 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. preCICE returns the
+// maximum admissible time-step size during the initialization.
+template <int dim, typename ParameterClass>
+double
+Adapter<dim, ParameterClass>::initialize(
+ const DoFHandler<dim> & dof_handler,
+ std::map<types::global_dof_index, double> &boundary_data,
+ const MappingQGeneric<dim> & mapping)
+{
+ Assert(dim > 1, ExcNotImplemented());
+ AssertDimension(dim, precice.getDimensions());
+
+ // In a first step, we get preCICE specific IDs from preCICE and store them in
+ // the respective variables. Later, they are used for data transfer.
+ mesh_id = precice.getMeshID(mesh_name);
+ read_data_id = precice.getDataID(read_data_name, mesh_id);
+
+
+ // 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()`
+ std::set<types::boundary_id> couplingBoundary;
+ couplingBoundary.insert(dealii_boundary_interface_id);
+
+ // The `ComponentMask()` might be important in case we deal with vector valued
+ // problems, because vector valued problems have a DoF for each component.
+ DoFTools::extract_boundary_dofs(dof_handler,
+ ComponentMask(),
+ coupling_dofs,
+ couplingBoundary);
+
+ // 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.
+ for (const auto i : coupling_dofs)
+ boundary_data[i] = 0.0;
+
+ // Since we deal with a scalar problem, the number of DoFs at the particular
+ // interface corresponds to the number of interface nodes.
+ n_interface_nodes = coupling_dofs.n_elements();
+
+ std::cout << "\t Number of coupling nodes: " << n_interface_nodes
+ << std::endl;
+
+ // 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.
+ std::vector<double> interface_nodes_positions;
+ interface_nodes_positions.reserve(dim * n_interface_nodes);
+
+ // 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.
+ read_data.resize(n_interface_nodes);
+ // The IDs are again filled by preCICE during the initializations.
+ interface_nodes_ids.resize(n_interface_nodes);
+
+ // The node location is obtained using `map_dofs_to_support_points()`.
+ std::map<types::global_dof_index, Point<dim>> support_points;
+ DoFTools::map_dofs_to_support_points(mapping, dof_handler, 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.
+ for (const auto element : coupling_dofs)
+ for (int i = 0; i < dim; ++i)
+ interface_nodes_positions.push_back(support_points[element][i]);
+
+ // Now we have all information to define the coupling mesh and pass the
+ // information to preCICE.
+ precice.setMeshVertices(mesh_id,
+ n_interface_nodes,
+ interface_nodes_positions.data(),
+ interface_nodes_ids.data());
+
+ // Then, we initialize preCICE internally calling the API function
+ // `initialize()`
+ const double max_delta_t = precice.initialize();
+
+
+ // read first coupling data from preCICE if available (i.e. deal.II is
+ // the second participant in a serial coupling scheme)
+ if (precice.isReadDataAvailable())
+ {
+ precice.readBlockScalarData(read_data_id,
+ n_interface_nodes,
+ interface_nodes_ids.data(),
+ read_data.data());
+
+ // After receiving the coupling data in `read_data`, we convert it to
+ // the std::map `boundary_data` which is later needed in order to apply
+ // Dirichlet boundary conditions
+ format_precice_to_dealii(boundary_data);
+ }
+
+ return max_delta_t;
+}
+
+
+// The function `advance()` is called in the main time loop after the
+// computation in each time step. Here,
+// coupling data is passed to and obtained from preCICE.
+template <int dim, typename ParameterClass>
+double
+Adapter<dim, ParameterClass>::advance(
+ std::map<types::global_dof_index, double> &boundary_data,
+ const double computed_timestep_length)
+{
+ // We specify the computed time-step length and pass it to preCICE. In
+ // return, preCICE tells us the maximum admissible time-step size our
+ // participant is allowed to compute in order to not exceed the next coupling
+ // time step.
+ const double max_delta_t = precice.advance(computed_timestep_length);
+
+ // As a next step, we obtain data, i.e. the boundary condition, from another
+ // participant. We have already all 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.
+ precice.readBlockScalarData(read_data_id,
+ n_interface_nodes,
+ interface_nodes_ids.data(),
+ read_data.data());
+
+ format_precice_to_dealii(boundary_data);
+
+ return max_delta_t;
+}
+
+
+
+// This function takes the std::vector obtained by preCICE in `read_data` 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` vector is exactly
+// the same as the order of the initially passed vertices coordinates.
+template <int dim, typename ParameterClass>
+void
+Adapter<dim, ParameterClass>::format_precice_to_dealii(
+ std::map<types::global_dof_index, double> &boundary_data) const
+{
+ // We already stored the coupling DoF indices in the `boundary_data` map, so
+ // that we can simply iterate over all keys in the map.
+ auto dof_component = boundary_data.begin();
+ for (int i = 0; i < n_interface_nodes; ++i)
+ {
+ AssertIndexRange(i, read_data.size());
+ boundary_data[dof_component->first] = read_data[i];
+ ++dof_component;
+ }
+}
+
+
+// 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.
+template <int dim>
+class CoupledLaplaceProblem
+{
+public:
+ CoupledLaplaceProblem();
+
+ void
+ run();
+
+private:
+ void
+ make_grid();
+ void
+ setup_system();
+ void
+ assemble_system();
+ void
+ solve();
+ void
+ output_results() const;
+
+ Triangulation<dim> triangulation;
+ FE_Q<dim> fe;
+ DoFHandler<dim> dof_handler;
+ MappingQ1<dim> mapping;
+
+ SparsityPattern sparsity_pattern;
+ SparseMatrix<double> system_matrix;
+
+ Vector<double> solution;
+ Vector<double> old_solution;
+ Vector<double> system_rhs;
+
+ // 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.
+ std::map<types::global_dof_index, double> boundary_data;
+ CouplingParamters parameters;
+ const types::boundary_id interface_boundary_id;
+ Adapter<dim, CouplingParamters> adapter;
+
+ // The time-step size delta_t is the acutual 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.
+ double delta_t;
+ double precice_delta_t;
+ const double solver_delta_t = 0.1;
+ unsigned int time_step = 0;
+};
+
+
+
+template <int dim>
+class RightHandSide : public Function<dim>
+{
+public:
+ virtual double
+ value(const Point<dim> &p, const unsigned int component = 0) const override;
+};
+
+
+
+template <int dim>
+class BoundaryValues : public Function<dim>
+{
+public:
+ virtual double
+ value(const Point<dim> &p, const unsigned int component = 0) const override;
+};
+
+template <int dim>
+double
+RightHandSide<dim>::value(const Point<dim> &p,
+ const unsigned int /*component*/) const
+{
+ double return_value = 0.0;
+ for (unsigned int i = 0; i < dim; ++i)
+ return_value += 4.0 * std::pow(p(i), 4.0);
+
+ return return_value;
+}
+
+
+template <int dim>
+double
+BoundaryValues<dim>::value(const Point<dim> &p,
+ const unsigned int /*component*/) const
+{
+ return p.square();
+}
+
+
+
+template <int dim>
+CoupledLaplaceProblem<dim>::CoupledLaplaceProblem()
+ : fe(1)
+ , dof_handler(triangulation)
+ , interface_boundary_id(1)
+ , adapter(parameters, interface_boundary_id)
+{}
+
+
+template <int dim>
+void
+CoupledLaplaceProblem<dim>::make_grid()
+{
+ GridGenerator::hyper_cube(triangulation, -1, 1);
+ triangulation.refine_global(4);
+
+ for (const auto &cell : triangulation.active_cell_iterators())
+ for (const auto face : cell->face_iterators())
+ {
+ // We choose the boundary in positive x direction for the
+ // interface coupling.
+ if (face->at_boundary() && (face->center()[0] == 1))
+ face->set_boundary_id(interface_boundary_id);
+ }
+
+ std::cout << " Number of active cells: " << triangulation.n_active_cells()
+ << std::endl
+ << " Total number of cells: " << triangulation.n_cells()
+ << std::endl;
+}
+
+
+template <int dim>
+void
+CoupledLaplaceProblem<dim>::setup_system()
+{
+ dof_handler.distribute_dofs(fe);
+
+ std::cout << " Number of degrees of freedom: " << dof_handler.n_dofs()
+ << std::endl;
+
+ DynamicSparsityPattern dsp(dof_handler.n_dofs());
+ DoFTools::make_sparsity_pattern(dof_handler, dsp);
+ sparsity_pattern.copy_from(dsp);
+
+ system_matrix.reinit(sparsity_pattern);
+
+ solution.reinit(dof_handler.n_dofs());
+ old_solution.reinit(dof_handler.n_dofs());
+ system_rhs.reinit(dof_handler.n_dofs());
+}
+
+
+
+template <int dim>
+void
+CoupledLaplaceProblem<dim>::assemble_system()
+{
+ // Reset global structures
+ system_rhs = 0;
+ system_matrix = 0;
+ // Update old solution values
+ old_solution = solution;
+
+ QGauss<dim> quadrature_formula(fe.degree + 1);
+
+ RightHandSide<dim> right_hand_side;
+
+ FEValues<dim> fe_values(fe,
+ quadrature_formula,
+ update_values | update_gradients |
+ update_quadrature_points | update_JxW_values);
+
+ const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
+
+ FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
+ Vector<double> cell_rhs(dofs_per_cell);
+ std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
+ // The solution values from previous time steps are stored for each quadrature
+ // point
+ std::vector<double> local_values_old_solution(fe_values.n_quadrature_points);
+
+ for (const auto &cell : dof_handler.active_cell_iterators())
+ {
+ fe_values.reinit(cell);
+ cell_matrix = 0;
+ cell_rhs = 0;
+ // Get the local values from the `fe_values' object
+ fe_values.get_function_values(old_solution, local_values_old_solution);
+
+ // The system matrix contains additionally a mass matrix due to the time
+ // discretization. The RHS has contributions from the old solution values.
+ for (const unsigned int q_index : fe_values.quadrature_point_indices())
+ for (const unsigned int i : fe_values.dof_indices())
+ {
+ for (const unsigned int j : fe_values.dof_indices())
+ cell_matrix(i, j) +=
+ ((fe_values.shape_value(i, q_index) * // phi_i(x_q)
+ fe_values.shape_value(j, q_index)) + // phi_j(x_q)
+ (delta_t * // delta t
+ fe_values.shape_grad(i, q_index) * // grad phi_i(x_q)
+ fe_values.shape_grad(j, q_index))) * // grad phi_j(x_q)
+ fe_values.JxW(q_index); // dx
+
+ const auto x_q = fe_values.quadrature_point(q_index);
+ const auto &local_value = local_values_old_solution[q_index];
+ cell_rhs(i) += ((delta_t * // delta t
+ fe_values.shape_value(i, q_index) * // phi_i(x_q)
+ right_hand_side.value(x_q)) + // f(x_q)
+ fe_values.shape_value(i, q_index) *
+ local_value) * // phi_i(x_q)*val
+ fe_values.JxW(q_index); // dx
+ }
+
+ // Copy local to global
+ cell->get_dof_indices(local_dof_indices);
+ for (const unsigned int i : fe_values.dof_indices())
+ {
+ for (const unsigned int j : fe_values.dof_indices())
+ system_matrix.add(local_dof_indices[i],
+ local_dof_indices[j],
+ cell_matrix(i, j));
+
+ system_rhs(local_dof_indices[i]) += cell_rhs(i);
+ }
+ }
+ {
+ // At first, we apply the Dirichlet boundary condition from step-4, as
+ // usual.
+ std::map<types::global_dof_index, double> boundary_values;
+ VectorTools::interpolate_boundary_values(dof_handler,
+ 0,
+ BoundaryValues<dim>(),
+ boundary_values);
+ MatrixTools::apply_boundary_values(boundary_values,
+ system_matrix,
+ solution,
+ system_rhs);
+ }
+ {
+ // Afterwards, we apply the coupling boundary condition. The `boundary_data`
+ // has already been filled by preCICE.
+ MatrixTools::apply_boundary_values(boundary_data,
+ system_matrix,
+ solution,
+ system_rhs);
+ }
+}
+
+
+
+template <int dim>
+void
+CoupledLaplaceProblem<dim>::solve()
+{
+ SolverControl solver_control(1000, 1e-12);
+ SolverCG<Vector<double>> solver(solver_control);
+ solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity());
+
+ std::cout << " " << solver_control.last_step()
+ << " CG iterations needed to obtain convergence." << std::endl;
+}
+
+
+
+template <int dim>
+void
+CoupledLaplaceProblem<dim>::output_results() const
+{
+ DataOut<dim> data_out;
+
+ data_out.attach_dof_handler(dof_handler);
+ data_out.add_data_vector(solution, "solution");
+
+ data_out.build_patches(mapping);
+
+ std::ofstream output("solution-" + std::to_string(time_step) + ".vtk");
+ data_out.write_vtk(output);
+}
+
+
+
+template <int dim>
+void
+CoupledLaplaceProblem<dim>::run()
+{
+ std::cout << "Solving problem in " << dim << " space dimensions."
+ << std::endl;
+
+ make_grid();
+ setup_system();
+
+ // After we set up the system, we initialize preCICE using the functionalities
+ // of the Adapter. preCICE returns the maximum admissible time-step size,
+ // which needs to be compared to our desired solver time-step size.
+ precice_delta_t = adapter.initialize(dof_handler, boundary_data, mapping);
+ delta_t = std::min(precice_delta_t, solver_delta_t);
+
+ // preCICE steers the coupled simulation: `isCouplingOngoing` is
+ // used to synchronize the end of the simulation with the coupling partner
+ while (adapter.precice.isCouplingOngoing())
+ {
+ // The time step number is solely used to generate unique output files
+ ++time_step;
+ // In the time loop, we assemble the coupled system and solve it as
+ // usual.
+ assemble_system();
+ solve();
+
+ // 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 and obtain new data from preCICE. Here, we simply
+ // obtain new data from preCICE, so from the other participant. As before,
+ // we obtain a maximum time-step size and compare it against the desired
+ // solver time-step size.
+ precice_delta_t = adapter.advance(boundary_data, delta_t);
+ delta_t = std::min(precice_delta_t, solver_delta_t);
+
+ // 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`.
+ if (adapter.precice.isTimeWindowComplete())
+ output_results();
+ }
+}
+
+
+
+int
+main()
+{
+ CoupledLaplaceProblem<2> laplace_problem;
+ laplace_problem.run();
+
+ return 0;
+}
projects / code-gallery.git / commitdiff
