23#ifndef INCLUDE_EXADG_FLUID_STRUCTURE_INTERACTION_ACCELERATION_SCHEMES_PARTITIONED_SOLVER_H_
24#define INCLUDE_EXADG_FLUID_STRUCTURE_INTERACTION_ACCELERATION_SCHEMES_PARTITIONED_SOLVER_H_
27#include <exadg/fluid_structure_interaction/acceleration_schemes/linear_algebra.h>
28#include <exadg/fluid_structure_interaction/acceleration_schemes/parameters.h>
29#include <exadg/fluid_structure_interaction/single_field_solvers/fluid.h>
30#include <exadg/fluid_structure_interaction/single_field_solvers/structure.h>
33#include <exadg/utilities/print_solver_results.h>
34#include <exadg/utilities/timer_tree.h>
40template<
int dim,
typename Number>
44 typedef dealii::LinearAlgebra::distributed::Vector<Number> VectorType;
54 solve(std::function<
void(VectorType &, VectorType
const &,
unsigned int)>
const &
55 apply_dirichlet_neumann_scheme);
58 print_iterations(dealii::ConditionalOStream
const & pcout)
const;
60 std::shared_ptr<TimerTree>
65 check_convergence(VectorType
const & residual)
const;
68 print_solver_info_header(
unsigned int const iteration)
const;
71 print_solver_info_converged(
unsigned int const iteration)
const;
76 dealii::ConditionalOStream pcout;
78 std::shared_ptr<SolverFluid<dim, Number>> fluid;
79 std::shared_ptr<SolverStructure<dim, Number>> structure;
82 std::vector<std::shared_ptr<std::vector<VectorType>>> D_history, R_history, Z_history;
85 std::shared_ptr<TimerTree> timer_tree;
91 std::pair<unsigned int, unsigned long long> partitioned_iterations;
94template<
int dim,
typename Number>
96 MPI_Comm
const & comm)
97 : parameters(parameters),
98 pcout(std::cout, dealii::Utilities::MPI::this_mpi_process(comm) == 0),
99 partitioned_iterations({0, 0})
101 timer_tree = std::make_shared<TimerTree>();
104template<
int dim,
typename Number>
106PartitionedSolver<dim, Number>::setup(std::shared_ptr<SolverFluid<dim, Number>> fluid_,
107 std::shared_ptr<SolverStructure<dim, Number>> structure_)
110 structure = structure_;
113template<
int dim,
typename Number>
115PartitionedSolver<dim, Number>::check_convergence(VectorType
const & residual)
const
117 double const residual_norm = residual.l2_norm();
118 double const ref_norm_abs = std::sqrt(structure->pde_operator->get_number_of_dofs());
119 double const ref_norm_rel = structure->time_integrator->get_velocity_np().l2_norm() *
120 structure->time_integrator->get_time_step_size();
122 bool const converged = (residual_norm < parameters.abs_tol * ref_norm_abs) or
123 (residual_norm < parameters.rel_tol * ref_norm_rel);
128template<
int dim,
typename Number>
130PartitionedSolver<dim, Number>::print_solver_info_header(
unsigned int const iteration)
const
132 if(fluid->time_integrator->print_solver_info())
135 <<
"======================================================================" << std::endl
136 <<
" Partitioned FSI: iteration counter = " << std::left << std::setw(8) << iteration
138 <<
"======================================================================" << std::endl;
142template<
int dim,
typename Number>
144PartitionedSolver<dim, Number>::print_solver_info_converged(
unsigned int const iteration)
const
146 if(fluid->time_integrator->print_solver_info())
149 <<
"Partitioned FSI iteration converged in " << iteration <<
" iterations." << std::endl;
153template<
int dim,
typename Number>
155PartitionedSolver<dim, Number>::print_iterations(dealii::ConditionalOStream
const & pcout)
const
157 std::vector<std::string> names;
158 std::vector<double> iterations_avg;
160 names = {
"Partitioned iterations"};
161 iterations_avg.resize(1);
163 (double)partitioned_iterations.second / std::max(1.0, (
double)partitioned_iterations.first);
165 print_list_of_iterations(pcout, names, iterations_avg);
168template<
int dim,
typename Number>
169std::shared_ptr<TimerTree>
170PartitionedSolver<dim, Number>::get_timings()
const
175template<
int dim,
typename Number>
177PartitionedSolver<dim, Number>::solve(
178 std::function<
void(VectorType &, VectorType
const &,
unsigned int)>
const &
179 apply_dirichlet_neumann_scheme)
185 if(parameters.acceleration_method == AccelerationMethod::Aitken)
188 structure->pde_operator->initialize_dof_vector(r_old);
189 structure->pde_operator->initialize_dof_vector(d);
191 bool converged =
false;
193 while(not(converged) and k < parameters.partitioned_iter_max)
195 print_solver_info_header(k);
198 structure->time_integrator->extrapolate_displacement_to_np(d);
200 d = structure->time_integrator->get_displacement_np();
202 VectorType d_tilde(d);
203 apply_dirichlet_neumann_scheme(d_tilde, d, k);
206 VectorType r = d_tilde;
208 converged = check_convergence(r);
218 omega = parameters.omega_init;
222 VectorType delta_r = r;
223 delta_r.add(-1.0, r_old);
224 omega *= -(r_old * delta_r) / delta_r.norm_sqr();
230 structure->time_integrator->set_displacement(d);
232 timer_tree->insert({
"Aitken"}, timer.wall_time());
239 else if(parameters.acceleration_method == AccelerationMethod::IQN_ILS)
241 std::shared_ptr<std::vector<VectorType>> D, R;
242 D = std::make_shared<std::vector<VectorType>>();
243 R = std::make_shared<std::vector<VectorType>>();
245 VectorType d, d_tilde, d_tilde_old, r, r_old;
246 structure->pde_operator->initialize_dof_vector(d);
247 structure->pde_operator->initialize_dof_vector(d_tilde);
248 structure->pde_operator->initialize_dof_vector(d_tilde_old);
249 structure->pde_operator->initialize_dof_vector(r);
250 structure->pde_operator->initialize_dof_vector(r_old);
252 unsigned int const q = parameters.reused_time_steps;
253 unsigned int const n = fluid->time_integrator->get_number_of_time_steps();
255 bool converged =
false;
256 while(not(converged) and k < parameters.partitioned_iter_max)
258 print_solver_info_header(k);
261 structure->time_integrator->extrapolate_displacement_to_np(d);
263 d = structure->time_integrator->get_displacement_np();
265 apply_dirichlet_neumann_scheme(d_tilde, d, k);
270 converged = check_convergence(r);
278 if(k == 0 and (q == 0 or n == 0))
280 d.add(parameters.omega_init, r);
287 VectorType delta_d_tilde = d_tilde;
288 delta_d_tilde.add(-1.0, d_tilde_old);
289 D->push_back(delta_d_tilde);
291 VectorType delta_r = r;
292 delta_r.add(-1.0, r_old);
293 R->push_back(delta_r);
297 std::vector<VectorType> Q = *R;
298 for(
auto R_q : R_history)
299 for(auto delta_r : *R_q)
300 Q.push_back(delta_r);
301 std::vector<VectorType> D_all = *D;
302 for(
auto D_q : D_history)
303 for(auto delta_d : *D_q)
304 D_all.push_back(delta_d);
306 AssertThrow(D_all.size() == Q.size(),
307 dealii::ExcMessage(
"D, Q vectors must have same size."));
309 unsigned int const k_all = Q.size();
313 Matrix<Number> U(k_all);
314 compute_QR_decomposition(Q, U);
316 std::vector<Number> rhs(k_all, 0.0);
317 for(
unsigned int i = 0; i < k_all; ++i)
318 rhs[i] = -Number(Q[i] * r);
321 std::vector<Number> alpha(k_all, 0.0);
322 backward_substitution(U, alpha, rhs);
326 for(
unsigned int i = 0; i < k_all; ++i)
327 d.add(alpha[i], D_all[i]);
331 d.add(parameters.omega_init, r);
335 d_tilde_old = d_tilde;
338 structure->time_integrator->set_displacement(d);
340 timer_tree->insert({
"IQN-ILS"}, timer.wall_time());
351 D_history.push_back(D);
352 R_history.push_back(R);
353 if(D_history.size() > q)
354 D_history.erase(D_history.begin());
355 if(R_history.size() > q)
356 R_history.erase(R_history.begin());
358 timer_tree->insert({
"IQN-ILS"}, timer.wall_time());
360 else if(parameters.acceleration_method == AccelerationMethod::IQN_IMVLS)
362 std::shared_ptr<std::vector<VectorType>> D, R;
363 D = std::make_shared<std::vector<VectorType>>();
364 R = std::make_shared<std::vector<VectorType>>();
366 std::vector<VectorType> B;
368 VectorType d, d_tilde, d_tilde_old, r, r_old, b, b_old;
369 structure->pde_operator->initialize_dof_vector(d);
370 structure->pde_operator->initialize_dof_vector(d_tilde);
371 structure->pde_operator->initialize_dof_vector(d_tilde_old);
372 structure->pde_operator->initialize_dof_vector(r);
373 structure->pde_operator->initialize_dof_vector(r_old);
374 structure->pde_operator->initialize_dof_vector(b);
375 structure->pde_operator->initialize_dof_vector(b_old);
377 std::shared_ptr<Matrix<Number>> U;
378 std::vector<VectorType> Q;
380 unsigned int const q = parameters.reused_time_steps;
381 unsigned int const n = fluid->time_integrator->get_number_of_time_steps();
383 bool converged =
false;
384 while(not(converged) and k < parameters.partitioned_iter_max)
386 print_solver_info_header(k);
389 structure->time_integrator->extrapolate_displacement_to_np(d);
391 d = structure->time_integrator->get_displacement_np();
393 apply_dirichlet_neumann_scheme(d_tilde, d, k);
398 converged = check_convergence(r);
407 inv_jacobian_times_residual(b, D_history, R_history, Z_history, r);
409 if(k == 0 and (q == 0 or n == 0))
411 d.add(parameters.omega_init, r);
421 VectorType delta_d_tilde = d_tilde;
422 delta_d_tilde.add(-1.0, d_tilde_old);
423 D->push_back(delta_d_tilde);
425 VectorType delta_r = r;
426 delta_r.add(-1.0, r_old);
427 R->push_back(delta_r);
429 VectorType delta_b = delta_d_tilde;
430 delta_b.add(1.0, b_old);
431 delta_b.add(-1.0, b);
432 B.push_back(delta_b);
435 U = std::make_shared<Matrix<Number>>(k);
437 compute_QR_decomposition(Q, *U);
439 std::vector<Number> rhs(k, 0.0);
440 for(
unsigned int i = 0; i < k; ++i)
441 rhs[i] = -Number(Q[i] * r);
444 std::vector<Number> alpha(k, 0.0);
445 backward_substitution(*U, alpha, rhs);
447 for(
unsigned int i = 0; i < k; ++i)
448 d.add(alpha[i], B[i]);
452 d_tilde_old = d_tilde;
456 structure->time_integrator->set_displacement(d);
458 timer_tree->insert({
"IQN-IMVLS"}, timer.wall_time());
469 D_history.push_back(D);
470 R_history.push_back(R);
471 if(D_history.size() > q)
472 D_history.erase(D_history.begin());
473 if(R_history.size() > q)
474 R_history.erase(R_history.begin());
477 std::shared_ptr<std::vector<VectorType>> Z;
478 Z = std::make_shared<std::vector<VectorType>>();
480 backward_substitution_multiple_rhs(*U, *Z, Q);
481 Z_history.push_back(Z);
482 if(Z_history.size() > q)
483 Z_history.erase(Z_history.begin());
485 timer_tree->insert({
"IQN-IMVLS"}, timer.wall_time());
489 AssertThrow(
false, dealii::ExcMessage(
"This AccelerationMethod is not implemented."));
492 partitioned_iterations.first += 1;
493 partitioned_iterations.second += k;
495 print_solver_info_converged(k);
Definition partitioned_solver.h:42
Definition structure.h:38
Definition parameters.h:38