coupled_laplace_problem.h

Go to the documentation of this file.

) 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())
 *         {
 * @endcode
 * 
 * We choose the boundary in positive x direction for the
 * interface coupling.
 * 
 * @code
 *           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()
 *   {
 * @endcode
 * 
 * Reset global structures
 * 
 * @code
 *     system_rhs    = 0;
 *     system_matrix = 0;
 * @endcode
 * 
 * Update old solution values
 * 
 * @code
 *     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);
 * @endcode
 * 
 * The solution values from previous time steps are stored for each quadrature
 * point
 * 
 * @code
 *     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;
 * @endcode
 * 
 * Get the local values from the `fe_values' object
 * 
 * @code
 *         fe_values.get_function_values(old_solution, local_values_old_solution);
 *   
 * @endcode
 * 
 * The system matrix contains additionally a mass matrix due to the time
 * discretization. The RHS has contributions from the old solution values.
 * 
 * @code
 *         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
 *             }
 *   
 * @endcode
 * 
 * Copy local to global
 * 
 * @code
 *         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);
 *           }
 *       }
 *     {
 * @endcode
 * 
 * At first, we apply the Dirichlet boundary condition from @ref step_4 "step-4", as
 * usual.
 * 
 * @code
 *       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);
 *     }
 *     {
 * @endcode
 * 
 * Afterwards, we apply the coupling boundary condition. The `boundary_data`
 * has already been filled by preCICE.
 * 
 * @code
 *       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();
 *   
 * @endcode
 * 
 * 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.
 * 
 * @code
 *     precice_delta_t = adapter.initialize(dof_handler, boundary_data, mapping);
 *     delta_t         = std::min(precice_delta_t, solver_delta_t);
 *   
 * @endcode
 * 
 * preCICE steers the coupled simulation: `isCouplingOngoing` is
 * used to synchronize the end of the simulation with the coupling partner
 * 
 * @code
 *     while (adapter.precice.isCouplingOngoing())
 *       {
 * @endcode
 * 
 * The time step number is solely used to generate unique output files
 * 
 * @code
 *         ++time_step;
 * @endcode
 * 
 * In the time loop, we assemble the coupled system and solve it as
 * usual.
 * 
 * @code
 *         assemble_system();
 *         solve();
 *   
 * @endcode
 * 
 * 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.
 * 
 * @code
 *         precice_delta_t = adapter.advance(boundary_data, delta_t);
 *         delta_t         = std::min(precice_delta_t, solver_delta_t);
 *   
 * @endcode
 * 
 * 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`.
 * 
 * @code
 *         if (adapter.precice.isTimeWindowComplete())
 *           output_results();
 *       }
 *   }
 *   
 *   
 *   
 *   int
 *   main()
 *   {
 *     CoupledLaplaceProblem<2> laplace_problem;
 *     laplace_problem.run();
 *   
 *     return 0;
 *   }
 * @endcode
 
 
<a name="ann-fancy_boundary_condition.cc"></a>
<h1>Annotated version of fancy_boundary_condition.cc</h1>
 * 
 * 
 * 
 * This program does not use any deal.II functionality and depends only on
 * preCICE and the standard libraries.
 * 
 * @code
 *   #include <precice/SolverInterface.hpp>
 *   
 *   #include <iostream>
 *   #include <sstream>
 *   
 * @endcode
 * 
 * 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
 * 
 * @code
 *   void
 *   define_boundary_values(std::vector<double> &boundary_data,
 *                          const double         time,
 *                          const double         end_time)
 *   {
 * @endcode
 * 
 * Scale the current time value
 * 
 * @code
 *     const double relative_time = time / end_time;
 * @endcode
 * 
 * Define the amplitude. Values run from -0.5 to 0.5
 * 
 * @code
 *     const double amplitude = (relative_time - 0.5);
 * @endcode
 * 
 * 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.
 * 
 * @code
 *     const double n_elements = boundary_data.size();
 *     const double right_zero = boundary_data.size() - 1;
 *     const double left_zero  = 0;
 *     const double offset     = 2;
 *     for (uint i = 0; i < n_elements; ++i)
 *       boundary_data[i] =
 *         -amplitude * ((i - left_zero) * (i - right_zero)) + offset;
 *   }
 *   
 *   
 *   int
 *   main()
 *   {
 *     std::cout << "Boundary participant: starting... \n";
 *   
 * @endcode
 * 
 * Configuration
 * 
 * @code
 *     const std::string configFileName("precice-config.xml");
 *     const std::string solverName("boundary-participant");
 *     const std::string meshName("boundary-mesh");
 *     const std::string dataWriteName("boundary-data");
 *   
 * @endcode
 * 
 * Adjust to MPI rank and size for parallel computation
 * 
 * @code
 *     const int commRank = 0;
 *     const int commSize = 1;
 *   
 *     precice::SolverInterface precice(solverName,
 *                                      configFileName,
 *                                      commRank,
 *                                      commSize);
 *   
 *     const int meshID           = precice.getMeshID(meshName);
 *     const int dimensions       = precice.getDimensions();
 *     const int numberOfVertices = 6;
 *   
 *     const int dataID = precice.getDataID(dataWriteName, meshID);
 *   
 * @endcode
 * 
 * Set up data structures
 * 
 * @code
 *     std::vector<double> writeData(numberOfVertices);
 *     std::vector<int>    vertexIDs(numberOfVertices);
 *     std::vector<double> vertices(numberOfVertices * dimensions);
 *   
 * @endcode
 * 
 * Define a boundary mesh
 * 
 * @code
 *     std::cout << "Boundary participant: defining boundary mesh \n";
 *     const double length = 2;
 *     const double xCoord = 1;
 *     const double deltaY = length / (numberOfVertices - 1);
 *     for (int i = 0; i < numberOfVertices; ++i)
 *       for (int j = 0; j < dimensions; ++j)
 *         {
 *           const unsigned int index = dimensions * i + j;
 * @endcode
 * 
 * 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.
 * 
 * @code
 *           if (j == 0)
 *             vertices[index] = xCoord;
 *           else
 *             vertices[index] = 1 - deltaY * i;
 *         }
 *   
 * @endcode
 * 
 * Pass the vertices to preCICE
 * 
 * @code
 *     precice.setMeshVertices(meshID,
 *                             numberOfVertices,
 *                             vertices.data(),
 *                             vertexIDs.data());
 *   
 * @endcode
 * 
 * initialize the Solverinterface
 * 
 * @code
 *     double dt = precice.initialize();
 *   
 * @endcode
 * 
 * Start time loop
 * 
 * @code
 *     const double end_time = 1;
 *     double       time     = 0;
 *     while (precice.isCouplingOngoing())
 *       {
 * @endcode
 * 
 * Generate new boundary data
 * 
 * @code
 *         define_boundary_values(writeData, time, end_time);
 *   
 *         {
 *           std::cout << "Boundary participant: writing coupling data \n";
 *           precice.writeBlockScalarData(dataID,
 *                                        numberOfVertices,
 *                                        vertexIDs.data(),
 *                                        writeData.data());
 *         }
 *   
 *         dt = precice.advance(dt);
 *         std::cout << "Boundary participant: advancing in time\n";
 *   
 *         time += dt;
 *       }
 *   
 *     std::cout << "Boundary participant: closing...\n";
 *   
 *     return 0;
 *   }
 * @endcode
 
 
*/