b2api
B2000++ API Reference Manual, VERSION 4.6
 
Loading...
Searching...
No Matches
b2dynamic_nonlinear_utile.H
1//------------------------------------------------------------------------
2// b2dynamic_nonlinear_utile.H --
3//
4//
5// written by Mathias Doreille
6// Neda Ebrahimi Pour <neda.ebrahimipour@dlr.de>
7//
8// (c) 2004-2012,2016 SMR Engineering & Development SA
9// 2502 Bienne, Switzerland
10//
11// (c) 2023 Deutsches Zentrum für Luft- und Raumfahrt (DLR) e.V.
12// Linder Höhe, 51147 Köln
13//
14// All Rights Reserved. Proprietary source code. The contents of
15// this file may not be disclosed to third parties, copied or
16// duplicated in any form, in whole or in part, without the prior
17// written permission of SMR.
18//------------------------------------------------------------------------
19
20#ifndef B2DYNAMIC_NONLINEAR_UTILE_H_
21#define B2DYNAMIC_NONLINEAR_UTILE_H_
22
23#include <cmath>
24
25#include "b2ppconfig.h"
26#include "model/b2boundary_condition.H"
27#include "model/b2case.H"
28#include "model/b2domain.H"
29#include "model/b2element.H"
30#include "model/b2model.H"
31#include "model/b2solution.H"
32#include "model/b2solver.H"
33#include "solvers/b2solver_utils.H"
34#include "solvers/b2static_nonlinear_utile.H"
35#include "utils/b2linear_algebra.H"
36#include "utils/b2logging.H"
37#include "utils/b2object.H"
38
39namespace {
40
41const int MAX_STEPS_DEFAULT = 99999;
42
43}
44
45namespace b2000 { namespace solver {
46
47template <typename T, typename MATRIX_FORMAT>
48class DynamicResidueFunction;
49
50template <typename T, typename MATRIX_FORMAT>
51class DynamicResidueFunction : public StaticResidueFunctionForTerminationTest<T> {
52public:
53 DynamicResidueFunction()
54 : logger(logging::get_logger("solver.static_nonlinear.ResidueFunction")) {}
55
56 virtual void init(Model& model, const AttributeList& attribute) = 0;
57
58 virtual void set_case(Case* case_) = 0;
59
60 virtual int get_order() const = 0;
61
62 virtual b2linalg::SparseMatrixConnectivityType get_dof_connectivity_type() const = 0;
63
64 virtual void save_solution(
65 b2linalg::Matrix<T, b2linalg::Mrectangle>& dof, const double time, const double dtime,
66 const RTable& attributes, bool end_of_stage = false,
67 bool unconverged_subcycle = false) = 0;
68
69 virtual void get_residue(
70 const b2linalg::Matrix<T, b2linalg::Mrectangle>& dof_red, const double time,
71 const double delta_time, const EquilibriumSolution equilibrium_solution,
72 b2linalg::Vector<T, b2linalg::Vdense>& residue,
73 const b2linalg::Vector<double, b2linalg::Vdense>& d_residue_d_dof_coeff,
74 b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>& d_residue_d_dofc,
75 b2linalg::Vector<T, b2linalg::Vdense>& d_residue_d_time, SolverHints* solver_hints,
76 CorrectionTerminationTest<T>* termination_test = 0) = 0;
77
78 virtual const b2linalg::Index& get_dof_field() = 0;
79
80 typedef ObjectTypeIncomplete<DynamicResidueFunction> type_t;
81 static type_t type;
82
83protected:
84 logging::Logger& logger;
85};
86
87template <typename T, typename MATRIX_FORMAT>
88typename DynamicResidueFunction<T, MATRIX_FORMAT>::type_t
89 DynamicResidueFunction<T, MATRIX_FORMAT>::type(
90 "DynamicResidueFunction", type_name<T, MATRIX_FORMAT>(),
91 StringList("DYNAMIC_RESIDUE_FUNCTION"), b2000::solver::module);
92
93class NordsieckParameter : public Object {
94public:
95 virtual void set_attribute(const AttributeList& attribute) {}
96
97 virtual void get_matrix(
98 int size, double h, int& order, b2linalg::Matrix<double, b2linalg::Mrectangle>& S,
99 b2linalg::Vector<double, b2linalg::Vdense>& d, double& truncated_error_coef,
100 double& stability_r_bound) = 0;
101
102 virtual int get_max_size() = 0;
103
104 virtual int get_min_size(int order) = 0;
105
106 virtual int get_order(int size) = 0;
107
108 typedef ObjectTypeIncomplete<NordsieckParameter> type_t;
109 static type_t type;
110
111protected:
112 static int factoriel(int i) {
113 if (i < 2) { return 1; }
114 int r = 1;
115 while (i != 1) { r *= i--; }
116 return r;
117 }
118
119 static double pow(int i, int j) {
120 if (j < 1) { return 1; }
121 double r = i;
122 while (--j != 0) { r *= i; }
123 return r;
124 }
125
126 static std::pair<int, double> get_order_and_truncated_error(
127 const b2linalg::Vector<double, b2linalg::Vdense_constref>& alpha,
128 const b2linalg::Vector<double, b2linalg::Vdense_constref>& beta) {
129 if (alpha.size() != beta.size()) { Exception() << THROW; }
130 double c = 1 - std::accumulate(&alpha[0], &alpha[0] + alpha.size(), 0.0);
131 if (std::abs(c) > 0.000001) { return std::pair<int, double>(-1, c); }
132 int order = 0;
133 while (size_t(order) <= alpha.size()) {
134 double cc = 0;
135 for (size_t i = 1; i != beta.size(); ++i) {
136 double i_pow = std::pow(double(i), order);
137 cc += i_pow * beta[i];
138 }
139 int q_fact = factoriel(order);
140 c = cc / q_fact;
141 cc = 0;
142 for (size_t i = 0; i != alpha.size(); ++i) {
143 double i_pow = std::pow(double(i + 1), order + 1);
144 cc += i_pow * alpha[i];
145 }
146 c -= cc / (q_fact * (order + 1));
147 if (order == 0) { c += beta[0]; }
148 if (order % 2 == 0) { c = -c; }
149 if (std::abs(c) > 0.000001) { return std::pair<int, double>(order, c); }
150 ++order;
151 }
152 return std::pair<int, double>(order, 0);
153 }
154};
155
156template <typename T, typename MATRIX_FORMAT>
157class NordsieckSolver : public Object {
158public:
159 enum NewtonMethod { conventional, modified, delayed_modified, modified_postponed };
160
161 NordsieckSolver()
162 : line_search(nullptr),
163 termination_test(0),
164 logger(logging::get_logger("solver.dynamic_nonlinear.nordsieck_solver")),
165 increment_number(0),
166 first_increment_number_of_stage(0),
167 max_increments(0),
168 max_time(0),
169 residue(0),
170 model(nullptr),
171 case_(nullptr),
172 nordsieck_parameter(nullptr),
173 explicit_nordsieck_parameter(nullptr),
174 h(0),
175 time(0),
176 nordsieck_size(0),
177 newton_method(modified_postponed),
178 newton_periodic_update(12),
179 newton_subperiod_residual_monitoring(12),
180 max_divergence(0),
181 new_integration_order(0),
182 last_step_with_line_search(false),
183 step_size_min(0),
184 step_size_max(0),
185 local_error_integration_tolerance_max_index(0),
186 max_newton_iterations(0),
187 corrector_nb_newton_iteration(0) {}
188
189 ~NordsieckSolver() override {
190 if (line_search) { line_search->free(); }
191 if (termination_test) { termination_test->free(); }
192 }
193
194 void set_attribute(Case& case__) {
195 const AttributeList& attribute = case__;
196 max_newton_iterations = attribute.get_int("MAX_NEWTON_ITERATIONS", 50);
197 max_divergence = attribute.get_int("MAX_DIVERGENCES", 4);
198 corrector_nb_newton_iteration = attribute.get_double("NB_ITERATION_CORRECTION", 0.5);
199 h = attribute.get_double("STEP_SIZE_INIT", 0.1);
200 step_size_min = attribute.get_double("STEP_SIZE_MIN", 1e-12);
201 step_size_max = attribute.get_double("STEP_SIZE_MAX", 1e100);
202 if (step_size_min < 0 || step_size_max < step_size_min) {
203 Exception() << "The minimum and/or maximum step size is "
204 "invalid."
205 << THROW;
206 }
207
208 if (h < step_size_min) {
209 Exception() << "The initial step size must be greater than or "
210 "equal to the minimum step size ("
211 << step_size_min << ")." << THROW;
212 }
213 if (h > step_size_max) {
214 Exception() << "The initial step size must be smaller than or "
215 "equal to the maximum step size ("
216 << step_size_max << ")." << THROW;
217 }
218 max_increments = size_t(case_->get_int("MAX_STEPS", MAX_STEPS_DEFAULT));
219
220 const std::pair<std::map<std::string, size_t>, b2linalg::Index>& df =
221 model->get_domain().get_dof_field();
222 const double default_leit = attribute.get_double("TOL_DYNAMIC", 1e-3);
223 local_error_integration_tolerance_max_index = 0;
224 for (std::map<std::string, size_t>::const_iterator i = df.first.begin();
225 i != df.first.end(); ++i) {
226 local_error_integration_tolerance[i->first] = std::pair<size_t, double>(
227 i->second, attribute.get_double("TOL_DYNAMIC." + i->first, default_leit));
228 local_error_integration_tolerance_max_index =
229 std::max(local_error_integration_tolerance_max_index, i->second);
230 }
231 ++local_error_integration_tolerance_max_index;
232 local_error_integration_tolerance["Lagrange_Multiplier"] = std::pair<size_t, double>(
233 local_error_integration_tolerance_max_index,
234 attribute.get_double("TOL_DYNAMIC.LAGRANGE_MULTIPLIER", 1e100));
235 ++local_error_integration_tolerance_max_index;
236
237 // Newton method.
238 {
239 std::string method_s = attribute.get_string("NEWTON_METHOD", "CONVENTIONAL");
240 if (method_s == "CONVENTIONAL") {
241 newton_method = conventional;
242 } else {
243 newton_periodic_update = attribute.get_int("NEWTON_PERIODIC_UPDATE", 100);
244 if (method_s == "MODIFIED") {
245 newton_method = modified;
246 } else if (method_s == "DELAYED_MODIFIED") {
247 newton_method = delayed_modified;
248 }
249 }
250 if (attribute.get_int("FULLNEWTON", 0) != 0) { newton_method = conventional; }
251 }
252
253 // Newton termination test
254 {
255 std::string termination_s =
256 attribute.get_string("CORRECTION_TERMINATION_TEST", "FLUX_NORMALISED");
257 if (termination_s != "") { termination_s += "_CORRECTION_TERMINATION_TEST"; }
258 termination_s = type_name<T>(termination_s);
259 typename CorrectionTerminationTest<T>::type_t* termination_t =
260 CorrectionTerminationTest<T>::type.get_subtype(termination_s, solver::module);
261 if (termination_test) { termination_test->free(); }
262 termination_test = termination_t->new_object(allocator);
263 termination_test->init(*model, attribute);
264 }
265
266 if (attribute.get_bool("LINE_SEARCH", false)) {
267 std::string line_search_s =
268 attribute.get_string("LINE_SEARCH_ALGORITHM", "STRONG_WOLFE");
269 if (line_search_s != "") { line_search_s += "_LINE_SEARCH"; }
270 typename LineSearch::type_t* line_search_t =
271 LineSearch::type.get_subtype(line_search_s, solver::module);
272 if (line_search) { line_search->free(); }
273 line_search = line_search_t->new_object(allocator);
274 line_search->init(*model, attribute);
275 } else {
276 line_search = nullptr;
277 }
278 }
279
280 void init(
281 Model& model_, DynamicResidueFunction<T, MATRIX_FORMAT>* residue_,
282 NordsieckParameter* implicit_nordsieck_parameter_,
283 NordsieckParameter* explicit_nordsieck_parameter_,
284 const b2linalg::Matrix<T, b2linalg::Mrectangle_constref>& initial_dof,
285 double initial_time, Case* case__) {
286 residue = residue_;
287 nordsieck_parameter = implicit_nordsieck_parameter_;
288 explicit_nordsieck_parameter = explicit_nordsieck_parameter_;
289 case_ = case__;
290 const int problem_order = residue->get_order();
291 if (problem_order < 1) {
292 Exception(
293 "The problem to solve is of order less than 1 in time - thus it is not dynamic.")
294 << THROW;
295 }
296 nordsieck_size = 2;
297 const size_t nb_dof = residue->get_number_of_dof();
298 Y.resize(nb_dof, problem_order * nordsieck_size);
299 Y.set_zero();
300 Y(b2linalg::Interval(0, initial_dof.size1()), b2linalg::Interval(0, initial_dof.size2())) =
301 initial_dof;
302 time = initial_time;
303 K.resize(nb_dof, nb_dof, residue_->get_dof_connectivity_type(), *case_);
304 increment_number = 0;
305 first_increment_number_of_stage = 0;
306 new_integration_order = true;
307 last_step_with_line_search = false;
308
309 model = &model_;
310 set_attribute(*case__);
311 }
312
313 bool next() {
314 Time start_analysis_time;
315 logger << logging::info << "Start increment " << increment_number << " of stage "
316 << case_->get_stage_id() << " for time " << time << " and time step " << h << "."
317 << logging::LOGFLUSH;
318
319 const std::string correction_termination_test_s =
320 case_->get_string("CORRECTION_TERMINATION_TEST", "NORMALISED_DOF_AND_RESIDUE");
321 bool converged = true;
322 double time_max = case_->get_stage_size();
323 bool is_refactorized_in_loop = false;
324 int max_nb_iter = 0;
325 const size_t number_of_dof = residue->get_number_of_dof();
326 const int problem_order = residue->get_order();
327 EquilibriumSolution equilibrium_solution(0, 1);
328 termination_test->new_step();
329 do {
330 // Get the Nordsieck parameters
331 b2linalg::Matrix<double, b2linalg::Mrectangle> explicit_S;
332 double explicit_truncated_error_coef;
333 int nordsieck_order;
334 b2linalg::Matrix<double, b2linalg::Mrectangle> S;
335 b2linalg::Vector<double, b2linalg::Vdense> d;
336 double truncated_error_coef;
337 double stability_r_bound;
338 explicit_nordsieck_parameter->get_matrix(
339 explicit_nordsieck_parameter->get_min_size(
340 nordsieck_parameter->get_order(nordsieck_size)),
341 h, nordsieck_order, explicit_S, d, explicit_truncated_error_coef,
342 stability_r_bound);
343 if (d[0] != 0) {
344 Exception() << "The explicit integration schema is not explicit." << THROW;
345 }
346 nordsieck_parameter->get_matrix(
347 nordsieck_size, h, nordsieck_order, S, d, truncated_error_coef,
348 stability_r_bound);
349
350 double milne_device_local_error_coef =
351 truncated_error_coef / (explicit_truncated_error_coef - truncated_error_coef);
352 double inv_d0 = 1.0 / d[0];
353 b2linalg::Vector<double, b2linalg::Vdense> d_residue_d_dof_coeff(problem_order + 1);
354 d_residue_d_dof_coeff[0] = 1;
355 d_residue_d_dof_coeff[1] = inv_d0;
356 for (int i = 1; i != problem_order; ++i) {
357 d_residue_d_dof_coeff[i + 1] = d_residue_d_dof_coeff[i] * inv_d0;
358 }
359
360 // Explicit part of implicit Nordsieck
361 b2linalg::Matrix<double, b2linalg::Mrectangle> explicit_dof(
362 number_of_dof, problem_order);
363 for (size_t j = 0; j != S.size2(); ++j) {
364 explicit_dof(b2linalg::Interval(0, problem_order)) +=
365 S(0, j) * Y(b2linalg::Interval(j * problem_order, (j + 1) * problem_order));
366 }
367
368 // Explicit schemas
369 b2linalg::Vector<T, b2linalg::Vdense> explicit_exp_dof(number_of_dof);
370 for (size_t i = 0; i != std::min(explicit_S.size2(), Y.size2()); ++i) {
371 explicit_exp_dof += explicit_S(0, i) * Y[i];
372 }
373
374 // Implicit Nordsieck corrector with MNR iteration
375 b2linalg::Vector<double, b2linalg::Vdense> res(number_of_dof);
376 b2linalg::Matrix<double, b2linalg::Mrectangle> dof(number_of_dof, problem_order + 1);
377 dof[0] = explicit_exp_dof;
378 if (logger.is_enabled_for(logging::data)) {
379 RTable attributes;
380 residue->save_solution(dof, time + h, h, attributes, false, true);
381 }
382 int nb_iter = 0;
383 b2linalg::Vector<T, b2linalg::Vindex1_constref> diag_K;
384 is_refactorized_in_loop = false;
385 b2linalg::Vector<T> delta_dof;
386 delta_dof = explicit_exp_dof - Y[0];
387 b2linalg::Vector<T> last_dof;
388 last_dof = Y[0];
389 SolverHints solver_hints;
390 typename CorrectionTerminationTest<T>::Termination result;
391 do {
392 termination_test->new_iteration();
393 termination_test->new_evaluation();
394
395 bool refactorization = (nb_iter == 0 && last_step_with_line_search)
396 || (newton_method == conventional)
397 || ((nb_iter + (newton_method == delayed_modified ? 0 : 1))
398 % newton_periodic_update
399 == 0)
400 || (!converged && nb_iter == 0) || new_integration_order;
401 new_integration_order = false;
402
403 for (int i = 0; i != problem_order; ++i) {
404 dof[i + 1] = inv_d0 * (dof[i] - explicit_dof[i]);
405 }
406
407 solver_hints.clear();
408 if (refactorization) {
409 residue->get_residue(
410 dof, time + h, h, equilibrium_solution, res, d_residue_d_dof_coeff, K,
411 b2linalg::Vector<T, b2linalg::Vdense>::null, &solver_hints,
412 termination_test->need_flux() ? termination_test : 0);
413 is_refactorized_in_loop = true;
414 if ((*solver_hints & SolverHints::diverge) != 0) {
415 logger << logging::info
416 << "elements/materials/boundary "
417 "conditions indicated divergence."
418 << logging::LOGFLUSH;
419 }
420 } else {
421 residue->get_residue(
422 dof, time + h, h, equilibrium_solution, res, d_residue_d_dof_coeff,
423 b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>::null,
424 b2linalg::Vector<T, b2linalg::Vdense>::null, &solver_hints,
425 termination_test->need_flux() ? termination_test : 0);
426 }
427 equilibrium_solution.input_level = 1;
428
429 if (logger.is_enabled_for(logging::data)) {
430 if (refactorization) {
431 logger << logging::data << "d_value_d_dof "
432 << logging::DBName(
433 "D_VALUE_D_DOF.GLOB", increment_number, nb_iter + 1,
434 case_->get_id())
435 << K << logging::LOGFLUSH;
436 }
437 logger << logging::data << "value"
438 << logging::DBName(
439 "VALUE.GLOB", increment_number, nb_iter + 1, case_->get_id())
440 << res << logging::LOGFLUSH;
441 }
442
443 if ((*solver_hints & SolverHints::terminate_analysis) != 0) {
444 return false;
445 } else if ((*solver_hints & SolverHints::diverge) != 0) {
446 result = CorrectionTerminationTest<T>::diverged;
447 } else if (nb_iter == 0) {
448 result = CorrectionTerminationTest<T>::unfinished;
449 } else {
450 result = termination_test->is_terminated(
451 nb_iter, res, delta_dof, 1.0, dof[0], last_dof, 0, 1, h, time + h,
452 *residue, K.get_diagonal(),
453 equilibrium_solution.distance_to_iterative_convergence, &solver_hints);
454 }
455 converged = result == CorrectionTerminationTest<T>::converged;
456
457 if (result == CorrectionTerminationTest<T>::unfinished) {
458 delta_dof = inverse(K) * res;
459 last_dof = dof[0];
460 dof[0] += -delta_dof;
461 }
462 ++nb_iter;
463 if (logger.is_enabled_for(logging::data)) {
464 RTable attributes;
465 residue->save_solution(dof, time + h, h, attributes, false, true);
466 }
467 } while (result == CorrectionTerminationTest<T>::unfinished);
468
469 max_nb_iter = std::max(nb_iter, max_nb_iter);
470
471 if (result == CorrectionTerminationTest<T>::diverged
472 || result == CorrectionTerminationTest<T>::max_iteration) {
473 if (is_refactorized_in_loop) {
474 h *= .5;
475 if (h < step_size_min) {
476 Exception() << "Cannot converge, the step size "
477 "is smaller than the minimum."
478 << THROW;
479 }
480 }
481 logger << logging::debug
482 << "The Newton iterations have diverged, new time step=" << h << "."
483 << logging::LOGFLUSH;
484 } else {
485 // Check the nordsieck local error
486 std::vector<T> local_error_est(local_error_integration_tolerance_max_index);
487 const b2linalg::Index& df_index = residue->get_dof_field();
488 for (size_t i = 0; i != explicit_exp_dof.size(); ++i) {
489 const double r = std::abs(dof(i, 0) - explicit_exp_dof[i]);
490 if (i < df_index.size()) {
491 local_error_est[df_index[i]] = std::max(local_error_est[df_index[i]], r);
492 } else {
493 local_error_est.back() = std::max(local_error_est.back(), r);
494 }
495 }
496 for (std::map<std::string, std::pair<size_t, double> >::const_iterator i =
497 local_error_integration_tolerance.begin();
498 i != local_error_integration_tolerance.end(); ++i) {
499 local_error_est[i->second.first] *= std::abs(milne_device_local_error_coef);
500 if (local_error_est[i->second.first] > i->second.second) {
501 // Rerun this step with a different step size
502 converged = false;
503 logger << logging::debug
504 << "Re-run this increment with a different step "
505 "size because the local time integration "
506 "estimation error "
507 << local_error_est[i->second.first]
508 << " exceeds the local error tolerance " << i->second.second
509 << " for the field " << i->first << "." << logging::LOGFLUSH;
510 }
511 }
512
513 if (converged) {
514 // Go to the next step
515 time += h;
516 for (int i = 0; i != problem_order; ++i) {
517 dof[i + 1] = inv_d0 * (dof[i] - explicit_dof[i]);
518 }
519 RTable attributes;
520 attributes.set("NUM_NEWTON_ITERATIONS", nb_iter);
521 Time end_analysis_time;
522 const double elapsed_time = double(end_analysis_time - start_analysis_time);
523 attributes.set("ELAPSED_TIME", elapsed_time);
524 residue->save_solution(dof, time, h, attributes, time >= time_max);
525 if (nordsieck_size < nordsieck_parameter->get_max_size()) {
526 ++nordsieck_size;
527 new_integration_order = true;
528 }
529 b2linalg::Matrix<T, b2linalg::Mrectangle> new_Y(
530 Y.size1(), problem_order * nordsieck_size);
531 new_Y(b2linalg::Interval(0, problem_order)) =
532 dof(b2linalg::Interval(0, problem_order));
533 for (size_t i = 1; i < S.size2(); ++i) {
534 new_Y(b2linalg::Interval(i * problem_order, (i + 1) * problem_order)) =
535 d[i] * dof(b2linalg::Interval(1, problem_order + 1));
536 for (size_t j = 0; j != S.size2(); ++j) {
537 new_Y(b2linalg::Interval(i * problem_order, (i + 1) * problem_order)) +=
538 S(i, j)
539 * Y(b2linalg::Interval(
540 j * problem_order, (j + 1) * problem_order));
541 }
542 }
543
544 Y.swap(new_Y);
545 termination_test->converged_step();
546 }
547
548 // Compute the new step size
549 double new_h_error = 1e100;
550 for (std::map<std::string, std::pair<size_t, double> >::const_iterator i =
551 local_error_integration_tolerance.begin();
552 i != local_error_integration_tolerance.end(); ++i) {
553 new_h_error = std::min(
554 new_h_error,
555 h
556 * std::pow(
557 i->second.second
558 / std::max(local_error_est[i->second.first] * 2, 1e-10),
559 1.0 / (nordsieck_order + 1)));
560 }
561 const double new_h_iter =
562 h
563 * std::pow(
564 corrector_nb_newton_iteration,
565 double(nb_iter - max_newton_iterations) / max_newton_iterations - 1);
566 const double bound_h_stability = h * stability_r_bound;
567 h = std::min(
568 step_size_max,
569 std::min(new_h_error, std::min(new_h_iter, bound_h_stability)));
570 h = std::max(step_size_min, h);
571 if (time + h * 1.00001 >= time_max) { h = time_max - time; }
572
573 if (h > new_h_error) {
574 Exception() << "The step size is smaller than the "
575 "minimum due to the local time error."
576 << THROW;
577 }
578 if (h > bound_h_stability) { Exception() << THROW; }
579
580 logger << logging::debug << "Computation of the new time step: error= (";
581 for (std::map<std::string, std::pair<size_t, double> >::const_iterator i =
582 local_error_integration_tolerance.begin();
583 i != local_error_integration_tolerance.end(); ++i) {
584 logger << i->first << ": " << local_error_est[i->second.first] << ", ";
585 }
586 logger << "), h_error= " << new_h_error << ", h_iter=" << new_h_iter
587 << ", h_stability=" << bound_h_stability << ", new time step=" << h << "."
588 << logging::LOGFLUSH;
589 }
590 equilibrium_solution.input_level = 0;
591 } while (!converged);
592
593 ++increment_number;
594 if (time >= time_max || (time_max - time) / time_max < 1e-10) {
595 if (!case_->next_stage()) { return false; }
596 time = 0;
597 residue->set_case(case_);
598 model->set_case(*case_);
599 set_attribute(*case_);
600 first_increment_number_of_stage = increment_number;
601 max_increments = size_t(case_->get_int("MAX_STEPS", MAX_STEPS_DEFAULT));
602 }
603
604 if (increment_number - first_increment_number_of_stage > max_increments) {
605 Exception() << "The number of increments exceeds the "
606 "maximum ("
607 << max_increments << ") for stage " << case_->get_stage_id() << "."
608 << THROW;
609 }
610
611 return true;
612 }
613
614 typedef ObjectTypeComplete<NordsieckSolver> type_t;
615 static type_t type;
616
617private:
618 struct Function : public LineSearch::Function {
619 Function(
620 DynamicResidueFunction<T, MATRIX_FORMAT>& residue_function_,
621 b2linalg::Vector<T, b2linalg::Vdense>& director_, b2linalg::Vector<double>& d_,
622 b2linalg::Matrix<double, b2linalg::Mrectangle>& d_test_,
623 b2linalg::Vector<T, b2linalg::Vdense>& r_test_, double time_, double time_h_,
624 double inv_d0_, b2linalg::Matrix<double, b2linalg::Mrectangle>& explicit_dof_)
625 : residue_function(residue_function_),
626 director(director_),
627 d(d_),
628 d_test(d_test_),
629 r_test(r_test_),
630 time(time_),
631 time_h(time_h_),
632 inv_d0(inv_d0_),
633 explicit_dof(explicit_dof_) {}
634
635 double operator()(double h) override {
636 d_test[0] = d - h * director;
637 for (size_t i = 0; i != explicit_dof.size2(); ++i) {
638 d_test[i + 1] = inv_d0 * (d_test[i] - explicit_dof[i]);
639 }
640 b2linalg::Vector<double, b2linalg::Vdense> d_residue_d_dof_coeff;
641 residue_function.get_residue(
642 d_test, time + time_h, time_h, EquilibriumSolution(1, -1), r_test,
643 d_residue_d_dof_coeff,
644 b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>::null,
645 b2linalg::Vector<T, b2linalg::Vdense>::null, 0);
646 double rh_norm = 0;
647 for (size_t i = 0; i != d_test.size1(); ++i) { rh_norm += norm(r_test[i] * d[i]); }
648 return rh_norm;
649 }
650
651 DynamicResidueFunction<T, MATRIX_FORMAT>& residue_function;
652 b2linalg::Vector<T, b2linalg::Vdense>& director;
653 b2linalg::Vector<double>& d;
654 b2linalg::Matrix<double, b2linalg::Mrectangle>& d_test;
655 b2linalg::Vector<T, b2linalg::Vdense>& r_test;
656 double time;
657 double time_h;
658 double inv_d0;
659 b2linalg::Matrix<double, b2linalg::Mrectangle>& explicit_dof;
660 };
661 LineSearch* line_search;
662 CorrectionTerminationTest<T>* termination_test;
663 Allocator allocator;
664
665 logging::Logger& logger;
666 size_t increment_number;
667 size_t first_increment_number_of_stage;
668 size_t max_increments;
669
670 // Problem definition
671 double max_time;
672 DynamicResidueFunction<T, MATRIX_FORMAT>* residue;
673 Model* model;
674 Case* case_;
675
676 // nordsieck data
677 NordsieckParameter* nordsieck_parameter;
678 NordsieckParameter* explicit_nordsieck_parameter;
679 b2linalg::Matrix<T, b2linalg::Mrectangle> Y;
680 double h;
681 double time;
682 int nordsieck_size;
683
684 // Jacobian matrix
685 b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible> K;
686 b2linalg::Vector<T, b2linalg::Vdense> v;
687
688 // Modified Newton parameter
689 NewtonMethod newton_method;
690 int newton_periodic_update;
691 int newton_subperiod_residual_monitoring;
692 int max_divergence;
693 bool new_integration_order;
694 bool last_step_with_line_search;
695
696 // Step size adaptation
697 double step_size_min;
698 double step_size_max;
699 size_t local_error_integration_tolerance_max_index;
700 std::map<std::string, std::pair<size_t, double> > local_error_integration_tolerance;
701 int max_newton_iterations;
702 double corrector_nb_newton_iteration;
703};
704
705template <typename T, typename MATRIX_FORMAT>
706typename NordsieckSolver<T, MATRIX_FORMAT>::type_t NordsieckSolver<T, MATRIX_FORMAT>::type(
707 "NordsieckSolver", type_name<T, MATRIX_FORMAT>(), StringList("NORDSIECK_SOLVER"),
708 b2000::solver::module);
709
710// =====================================================
711// Implementation
712// =====================================================
713
714class NordsieckBDFParameter : public NordsieckParameter {
715public:
716 virtual void get_BDF_coef(
717 b2linalg::Matrix<double, b2linalg::Mrectangle>& alpha,
718 b2linalg::Vector<double, b2linalg::Vdense>& beta,
719 b2linalg::Vector<double, b2linalg::Vdense>& stability_r_bound) const = 0;
720
721 int get_max_size() override {
722 if (list_precomputed_S.empty()) { init(); }
723 return list_precomputed_S.back().size1();
724 };
725
726 int get_min_size(int order) override {
727 if (list_precomputed_S.empty()) { init(); }
728 std::vector<int>::iterator i =
729 std::lower_bound(list_orders.begin(), list_orders.end(), order);
730 if (i == list_orders.end() || *i != order) { Exception() << THROW; }
731 return int(i - list_orders.begin()) + 1;
732 }
733
734 int get_order(int size) override {
735 if (list_precomputed_S.empty()) { init(); }
736 size -= 2;
737 if (size < 0 && size_t(size) >= list_precomputed_S.size()) { Exception() << THROW; }
738 return list_orders[size];
739 }
740
741 void get_matrix(
742 int size, double h, int& order, b2linalg::Matrix<double, b2linalg::Mrectangle>& S,
743 b2linalg::Vector<double, b2linalg::Vdense>& d, double& truncated_error_constant,
744 double& stability_r_bound) override {
745 if (list_precomputed_S.empty()) { init(); }
746 size -= 2;
747 if (size < 0 && size_t(size) >= list_precomputed_S.size()) { Exception() << THROW; }
748 S = list_precomputed_S[size];
749 d = h * list_d[size];
750 double hi = 1;
751 for (size_t i = 0; i < S.size1(); ++i) {
752 double hj = hi;
753 for (size_t j = 0; j < S.size2(); ++j) {
754 S(i, j) *= hj;
755 hj *= h;
756 }
757 d[i] *= hi;
758 hi /= h;
759 }
760
761 order = list_orders[size];
762 truncated_error_constant = list_truncated_error_constant[size];
763 stability_r_bound = list_stability_r_bound[size];
764 }
765
766 typedef ObjectTypeIncomplete<NordsieckBDFParameter, NordsieckParameter::type_t> type_t;
767 static type_t type;
768
769private:
770 void init();
771 std::vector<b2linalg::Matrix<double, b2linalg::Mrectangle> > list_precomputed_S;
772 std::vector<b2linalg::Vector<double, b2linalg::Vdense> > list_d;
773 std::vector<int> list_orders;
774 std::vector<double> list_truncated_error_constant;
775 std::vector<double> list_stability_r_bound;
776};
777
778class StandardBDFParameter : public NordsieckBDFParameter {
779public:
780 void set_attribute(const AttributeList& attribute) override {
781 max_size = attribute.get_int(
782 "MULTISTEP_INTEGRATION_ORDER", attribute.get_int("DYNAMIC_SOLVER_MAX_ORDER", 2));
783 }
784
785 void get_BDF_coef(
786 b2linalg::Matrix<double, b2linalg::Mrectangle>& alpha,
787 b2linalg::Vector<double, b2linalg::Vdense>& beta,
788 b2linalg::Vector<double, b2linalg::Vdense>& stability_r_bound) const override;
789
790 typedef ObjectTypeComplete<StandardBDFParameter, NordsieckBDFParameter::type_t> type_t;
791 static type_t type;
792
793private:
794 int max_size;
795};
796
797class NordsieckAdamsParameter : public NordsieckParameter {
798public:
799 virtual void get_adams_coef(
800 b2linalg::Matrix<double, b2linalg::Mrectangle>& beta,
801 b2linalg::Vector<double, b2linalg::Vdense>& stability_r_bound) const = 0;
802
803 int get_max_size() override {
804 if (list_precomputed_S.empty()) { init(); }
805 return list_precomputed_S.back().size1();
806 };
807
808 int get_min_size(int order) override {
809 if (list_precomputed_S.empty()) { init(); }
810 std::vector<int>::iterator i =
811 std::lower_bound(list_orders.begin(), list_orders.end(), order);
812 if (i == list_orders.end() || *i != order) { Exception() << THROW; }
813 return int(i - list_orders.begin()) + 1;
814 }
815
816 int get_order(int size) override {
817 if (list_precomputed_S.empty()) { init(); }
818 --size;
819 if (size < 0 && size_t(size) >= list_precomputed_S.size()) { Exception() << THROW; }
820 return list_orders[size];
821 }
822
823 void get_matrix(
824 int size, double h, int& order, b2linalg::Matrix<double, b2linalg::Mrectangle>& S,
825 b2linalg::Vector<double, b2linalg::Vdense>& d, double& truncated_error_constant,
826 double& stability_r_bound) override {
827 if (list_precomputed_S.empty()) { init(); }
828 --size;
829 if (size < 0 && size_t(size) >= list_precomputed_S.size()) { Exception() << THROW; }
830 S = list_precomputed_S[size];
831 d = h * list_d[size];
832 double hi = 1;
833 for (size_t i = 0; i < S.size1(); hi *= ++i / h) {
834 double hj = hi;
835 for (size_t j = 0; j < S.size2(); hj *= h / ++j) { S(i, j) *= hj; }
836 d[i] *= hi;
837 }
838
839 order = list_orders[size];
840 truncated_error_constant = list_truncated_error_constant[size];
841 stability_r_bound = list_stability_r_bound[size];
842 }
843
844 typedef ObjectTypeIncomplete<NordsieckAdamsParameter, NordsieckParameter::type_t> type_t;
845 static type_t type;
846
847private:
848 void init();
849 std::vector<b2linalg::Matrix<double, b2linalg::Mrectangle> > list_precomputed_S;
850 std::vector<b2linalg::Vector<double, b2linalg::Vdense> > list_d;
851 std::vector<int> list_orders;
852 std::vector<double> list_truncated_error_constant;
853 std::vector<double> list_stability_r_bound;
854};
855
856class AdamsBashforthParameter : public NordsieckAdamsParameter {
857public:
858 void set_attribute(const AttributeList& attribute) override {
859 max_size = attribute.get_int("DYNAMIC_SOLVER_MAX_ORDER", 4) + 2;
860 }
861
862 void get_adams_coef(
863 b2linalg::Matrix<double, b2linalg::Mrectangle>& beta,
864 b2linalg::Vector<double, b2linalg::Vdense>& stability_r_bound) const override;
865
866 typedef ObjectTypeComplete<AdamsBashforthParameter, NordsieckAdamsParameter::type_t> type_t;
867 static type_t type;
868
869private:
870 int max_size;
871};
872
873template <typename T, typename MATRIX_FORMAT = b2linalg::Msymmetric>
874class DynamicOrderNResidueFunction : public DynamicResidueFunction<T, MATRIX_FORMAT> {
875public:
876 typedef TypedNaturalBoundaryCondition<T, typename MATRIX_FORMAT::sparse_inversible> nbc_t;
877 typedef TypedEssentialBoundaryCondition<T> ebc_t;
878 typedef TypedModelReductionBoundaryCondition<T> mrbc_t;
879
880 DynamicOrderNResidueFunction()
881 : logger(logging::get_logger("solver.dynamic_nonlinear.residue")),
882 number_of_dof(0),
883 number_of_dof_red(0),
884 model(nullptr),
885 domain(nullptr),
886 mrbc(0),
887 ebc(0),
888 nbc(0),
889 solution(0),
890 constraint_has_penalty(false),
891 penalty_factor(0),
892 lagrange_mult_scale(1),
893 rcfo_only_spc(false) {}
894
895 void set_case(Case* case_) {}
896
897 void init(Model& model_, const AttributeList& attribute) {
898 model = &model_;
899 domain = &model->get_domain();
900 ebc = &dynamic_cast<ebc_t&>(model->get_essential_boundary_condition(ebc_t::type.get_subtype(
901 "TypedEssentialBoundaryConditionListComponent" + type_name<T>())));
902 mrbc = &dynamic_cast<mrbc_t&>(
903 model->get_model_reduction_boundary_condition(mrbc_t::type.get_subtype(
904 "TypedModelReductionBoundaryConditionListComponent" + type_name<T>())));
905 nbc = &dynamic_cast<nbc_t&>(model->get_natural_boundary_condition(nbc_t::type.get_subtype(
906 "TypedNaturalBoundaryConditionListComponent"
907 + type_name<T, typename MATRIX_FORMAT::sparse_inversible>())));
908 solution = &dynamic_cast<SolutionDynamicNonlinear<T>&>(model->get_solution());
909 number_of_dof = domain->get_number_of_dof();
910 nbc_sol.resize(number_of_dof);
911 fi.resize(number_of_dof);
912
913 constraint_has_penalty =
914 attribute.get_string("CONSTRAINT_METHOD", "REDUCTION") == "PENALTY";
915
916 constraint_has_penalty |=
917 attribute.get_string("NONLINEAR_CONSTRAINT_METHOD", "OTHER") == "PENALTY";
918
919 number_of_dof_red = mrbc->get_size(false, b2linalg::Vector<T, b2linalg::Vdense>::null, 1)
920 + (constraint_has_penalty ? 0 : ebc->get_size(false));
921
922 if (constraint_has_penalty) {
923 penalty_factor = attribute.get_double("CONSTRAINT_PENALTY", 1e10);
924 } else {
925 if (attribute.get_bool("REMOVE_DEPENDENT_CONSTRAINTS", false)) {
926 b2linalg::Matrix<T, b2linalg::Mcompressed_col> trans_d_l_d_dof;
927 ebc->get_nonlinear_value(
928 fi, 0, false, b2linalg::Vector<T, b2linalg::Vdense_ref>::null,
929 trans_d_l_d_dof, b2linalg::Vector<T, b2linalg::Vdense_ref>::null, 0);
930 trans_d_l_d_dof.get_dep_columns(ebc_to_remove, 3e-13, 1.5);
931 number_of_dof_red -= ebc_to_remove.size();
932 if (!ebc_to_remove.empty()) {
933 DynamicResidueFunction<T, MATRIX_FORMAT>::logger
934 << logging::info << "Remove " << ebc_to_remove.size()
935 << " dependent constraints." << logging::LOGFLUSH;
936 }
937 }
938 penalty_factor = attribute.get_double("CONSTRAINT_PENALTY", 1);
939 if (attribute.get_string("NONLINEAR_CONSTRAINT_METHOD", "OTHER") == "LAGRANGE") {
940 penalty_factor = 0;
941 }
942 lagrange_mult_scale = attribute.get_double("LAGRANGE_MULT_SCALE_PENALTY", 1.0);
943 }
944
945 rcfo_only_spc = attribute.get_bool("RCFO_ONLY_SPC", false);
946 if (rcfo_only_spc
947 && attribute.get_string("CONSTRAINT_METHOD", "REDUCTION") == "REDUCTION") {
948 Exception() << "It is not possible to compute the reaction "
949 "forces with rcfo_only_spc=yes and "
950 "constraint_method=reduction (which is the default). "
951 "Instead, set constraint_method either to lagrange, "
952 "augmented_lagrange, or penalty."
953 << THROW;
954 }
955 }
956
957 void mrbc_get_nonlinear_value(
958 const b2linalg::Vector<T, b2linalg::Vdense_constref>& dof, double time,
959 EquilibriumSolution equilibrium_solution, b2linalg::Vector<T, b2linalg::Vdense_ref> value,
960 b2linalg::Matrix<T, b2linalg::Mcompressed_col>& d_value_d_dof_trans,
961 b2linalg::Vector<T, b2linalg::Vdense_ref> d_value_d_time, SolverHints* solver_hints) {
962 mrbc->get_nonlinear_value(
963 dof, time, equilibrium_solution, value, d_value_d_dof_trans,
964 b2linalg::Matrix<T, b2linalg::Mrectangle_ref>(d_value_d_time), solver_hints);
965 }
966
967 size_t get_number_of_dof() {
968 return mrbc->get_size(false, b2linalg::Vector<T, b2linalg::Vdense>::null, 1)
969 + (constraint_has_penalty ? 0 : ebc->get_size(false) - ebc_to_remove.size());
970 }
971
972 size_t get_number_of_non_lagrange_dof() {
973 return mrbc->get_size(false, b2linalg::Vector<T, b2linalg::Vdense>::null, 1);
974 }
975
976 b2linalg::SparseMatrixConnectivityType get_dof_connectivity_type() const {
977 return domain->get_dof_connectivity_type();
978 }
979
980 int get_order() const { return domain->get_value_info().size() - 1; }
981
982 void save_solution(
983 b2linalg::Matrix<T, b2linalg::Mrectangle>& dof_red, const double time, const double dtime,
984 const RTable& attributes, bool end_of_stage = false, bool unconverged_subcycle = false) {
985 const size_t mrbc_size =
986 mrbc->get_size(false, b2linalg::Vector<T, b2linalg::Vdense>::null, 1);
987 b2linalg::Interval disp_range(0, mrbc_size);
988 b2linalg::Matrix<T, b2linalg::Mrectangle> dof(number_of_dof, dof_red.size2());
989 b2linalg::Matrix<T, b2linalg::Mcompressed_col> trans_d_r_d_dof;
990 b2linalg::Matrix<T, b2linalg::Mrectangle> d_r_d_time(number_of_dof, dof_red.size2() - 1);
991 mrbc->get_nonlinear_value(
992 dof_red[0](disp_range), time, true, dof[0], trans_d_r_d_dof, d_r_d_time, 0);
993 for (size_t j = 1; j < dof.size2(); ++j) {
994 dof[j] = transposed(trans_d_r_d_dof) * dof_red[j](disp_range);
995 dof[j] += d_r_d_time[j - 1];
996 }
997 if (rcfo_only_spc) {
998 // recompute fi using only the spc contribution
999
1000 b2linalg::Matrix<T, b2linalg::Mcompressed_col> trans_d_l_d_dof;
1001 b2linalg::Vector<T> l(ebc->get_size(false));
1002 ebc->get_nonlinear_value(
1003 dof[0], time, false, l, trans_d_l_d_dof,
1004 b2linalg::Vector<T, b2linalg::Vdense_ref>::null, 0);
1005 if (!ebc_to_remove.empty()) {
1006 l.remove(ebc_to_remove);
1007 trans_d_l_d_dof.remove_col(ebc_to_remove);
1008 }
1009
1010 if (constraint_has_penalty) {
1011 l *= -penalty_factor;
1012 } else {
1013 b2linalg::Interval lagrange_mult_range(mrbc_size, dof_red.size1());
1014 if (MATRIX_FORMAT::symmetric && penalty_factor != 0) {
1015 // We apply penalty method with a small factor to obtain a definite positive
1016 // matrix.
1017 l *= -penalty_factor;
1018 l += -lagrange_mult_scale * dof_red[0](lagrange_mult_range);
1019 } else {
1020 l = -lagrange_mult_scale * dof_red[0](lagrange_mult_range);
1021 }
1022 }
1023 for (size_t j = 0; j != trans_d_l_d_dof.size2(); ++j) {
1024 if (trans_d_l_d_dof.get_nb_nonzero(j) != 1) { l[j] = 0; }
1025 }
1026 fi = trans_d_l_d_dof * l;
1027 }
1028 if (unconverged_subcycle) {
1029 solution->set_solution(dof, time, nbc_sol, fi, unconverged_subcycle);
1030 return;
1031 }
1032 solution->new_cycle(end_of_stage);
1033 solution->set_solution(dof, time, nbc_sol, fi);
1034
1035 typename SolutionDynamicNonlinear<T>::gradient_container gc =
1036 solution->get_gradient_container();
1037 domain->set_dof(dof);
1038 EquilibriumSolution es(1, 0);
1039 solver_utile.get_elements_values(
1040 dof, time, dtime, es, false, trans_d_r_d_dof, d_r_d_time[0],
1041 b2linalg::Vector<double, b2linalg::Vdense>::null, *model,
1042 b2linalg::Vector<T, b2linalg::Vdense>::null,
1043 b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>::null,
1044 b2linalg::Vector<T, b2linalg::Vdense>::null, gc.get(), 0, logger, 0);
1045 std::string domain_state_id = solution->get_domain_state_id();
1046 if (!domain_state_id.empty()) { domain->save_state(domain_state_id); }
1047 solution->terminate_cycle(time, dtime, attributes);
1048 if (end_of_stage) { solution->terminate_stage(true); }
1049 }
1050
1051 void get_residue(
1052 const b2linalg::Matrix<T, b2linalg::Mrectangle>& dof_red, const double time,
1053 const double delta_time, const EquilibriumSolution equilibrium_solution,
1054 b2linalg::Vector<T, b2linalg::Vdense>& residue,
1055 const b2linalg::Vector<double, b2linalg::Vdense>& d_residue_d_dof_coeff,
1056 b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>& d_residue_d_dofc,
1057 b2linalg::Vector<T, b2linalg::Vdense>& d_residue_d_time, SolverHints* solver_hints,
1058 CorrectionTerminationTest<T>* termination_test = 0) {
1059 b2linalg::Interval disp_range(
1060 0, mrbc->get_size(false, b2linalg::Vector<T, b2linalg::Vdense>::null, 1));
1061 b2linalg::Interval lagrange_mult_range(
1062 mrbc->get_size(false, b2linalg::Vector<T, b2linalg::Vdense>::null, 1),
1063 number_of_dof_red);
1064 if (!residue.is_null()) { residue.set_zero(); }
1065 if (!d_residue_d_dofc.is_null()) { d_residue_d_dofc.set_zero_same_struct(); }
1066
1067 b2linalg::Matrix<T, b2linalg::Mcompressed_col> trans_d_r_d_dof;
1068 b2linalg::Matrix<T, b2linalg::Mrectangle> d_r_d_time(number_of_dof, dof_red.size2() - 1);
1069 b2linalg::Matrix<T, b2linalg::Mrectangle> dof(number_of_dof, dof_red.size2());
1070 mrbc->get_nonlinear_value(
1071 dof_red[0](disp_range), time, equilibrium_solution, dof[0], trans_d_r_d_dof,
1072 d_r_d_time, solver_hints);
1073 for (size_t j = 1; j < dof.size2(); ++j) {
1074 dof[j] = transposed(trans_d_r_d_dof) * dof_red[j](disp_range);
1075 dof[j] += d_r_d_time[j - 1];
1076 }
1077 b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible> d_fe_d_dof;
1078 b2linalg::Vector<T, b2linalg::Vdense> d_fe_d_time(number_of_dof);
1079 nbc_sol.set_zero();
1080 nbc->get_nonlinear_value(
1081 dof[0], time, equilibrium_solution,
1082 residue.is_null() ? b2linalg::Vector<T, b2linalg::Vdense_ref>::null
1083 : b2linalg::Vector<T, b2linalg::Vdense_ref>(nbc_sol),
1084 d_residue_d_dofc.is_null()
1085 ? b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>::null
1086 : d_fe_d_dof,
1087 d_residue_d_time.is_null() ? b2linalg::Vector<T, b2linalg::Vdense_ref>::null
1088 : b2linalg::Vector<T, b2linalg::Vdense_ref>(d_fe_d_time),
1089 solver_hints);
1090 if (termination_test && !residue.is_null()) { termination_test->add_flux(nbc_sol); }
1091 fi = -nbc_sol;
1092 d_fe_d_time = -d_fe_d_time;
1093 if (d_fe_d_dof.size1() != 0) {
1094 d_residue_d_dofc += -(trans_d_r_d_dof * d_fe_d_dof * transposed(trans_d_r_d_dof));
1095 }
1096
1097 b2linalg::Matrix<T, b2linalg::Mcompressed_col> trans_d_l_d_dof;
1098 b2linalg::Vector<T> l;
1099
1100 if (constraint_has_penalty) {
1101 if (!residue.is_null()) { l.resize(ebc->get_size(false)); }
1102 ebc->get_nonlinear_value(
1103 dof[0], time, equilibrium_solution,
1104 (residue.is_null() ? b2linalg::Vector<T, b2linalg::Vdense_ref>::null
1105 : b2linalg::Vector<T, b2linalg::Vdense_ref>(l)),
1106 trans_d_l_d_dof, b2linalg::Vector<T, b2linalg::Vdense_ref>::null, solver_hints);
1107 } else {
1108 if (ebc_to_remove.empty()) {
1109 ebc->get_nonlinear_value(
1110 dof[0], time, equilibrium_solution,
1111 (residue.is_null() ? b2linalg::Vector<T, b2linalg::Vdense_ref>::null
1112 : residue(lagrange_mult_range)),
1113 trans_d_l_d_dof,
1114 (d_residue_d_time.is_null() ? b2linalg::Vector<T, b2linalg::Vdense_ref>::null
1115 : d_residue_d_time(lagrange_mult_range)),
1116 solver_hints);
1117 } else {
1118 b2linalg::Vector<T> l(residue.is_null() ? 0 : ebc->get_size(false));
1119 ebc->get_nonlinear_value(
1120 dof[0], time, equilibrium_solution,
1121 (residue.is_null() ? b2linalg::Vector<T, b2linalg::Vdense_ref>::null
1122 : b2linalg::Vector<T, b2linalg::Vdense_ref>(l)),
1123 trans_d_l_d_dof,
1124 (d_residue_d_time.is_null() ? b2linalg::Vector<T, b2linalg::Vdense_ref>::null
1125 : d_residue_d_time(lagrange_mult_range)),
1126 solver_hints);
1127 if (!residue.is_null()) {
1128 l.remove(ebc_to_remove);
1129 residue(lagrange_mult_range) = l;
1130 }
1131 trans_d_l_d_dof.remove_col(ebc_to_remove);
1132 }
1133 if (!d_residue_d_time.is_null()) {
1134 d_residue_d_time(lagrange_mult_range) +=
1135 transposed(trans_d_l_d_dof) * d_r_d_time[0];
1136 }
1137 }
1138
1139 if (!d_residue_d_time.is_null() && d_fe_d_dof.size1()) {
1140 d_fe_d_time -= d_fe_d_dof * d_r_d_time[0];
1141 }
1142
1143 solver_utile.get_elements_values(
1144 dof, time, delta_time, equilibrium_solution, false, trans_d_r_d_dof, d_r_d_time[0],
1145 d_residue_d_dof_coeff, *model, fi, d_residue_d_dofc,
1146 d_residue_d_time.is_null() ? b2linalg::Vector<T, b2linalg::Vdense_ref>::null
1147 : b2linalg::Vector<T, b2linalg::Vdense_ref>(d_fe_d_time),
1148 0, solver_hints, logger, termination_test);
1149
1150 if (!residue.is_null()) {
1151 b2linalg::Vector<T> tmp = fi;
1152 if (constraint_has_penalty) {
1153 tmp += (trans_d_l_d_dof * l) * penalty_factor;
1154 } else {
1155 if (MATRIX_FORMAT::symmetric && penalty_factor != 0) {
1156 // We apply penalty method with a small factor to obtain a definite positive
1157 // matrix.
1158 b2linalg::Vector<T> tmp1;
1159 tmp1 = residue(lagrange_mult_range) * penalty_factor;
1160 tmp1 += dof_red[0](lagrange_mult_range) * lagrange_mult_scale;
1161 tmp += trans_d_l_d_dof * tmp1;
1162 } else {
1163 tmp +=
1164 lagrange_mult_scale * (trans_d_l_d_dof * dof_red[0](lagrange_mult_range));
1165 }
1166 residue(lagrange_mult_range) *= lagrange_mult_scale;
1167 }
1168 residue(disp_range) = trans_d_r_d_dof * tmp;
1169 }
1170 if (!d_residue_d_dofc.is_null()) {
1171 b2linalg::Matrix<T, b2linalg::Mcompressed_col> trans_LR;
1172 trans_LR = trans_d_r_d_dof * trans_d_l_d_dof;
1173 if (penalty_factor != 0 && (constraint_has_penalty || MATRIX_FORMAT::symmetric)) {
1174 // We apply penalty method but with a small factor to
1175 // obtain a definite positive matrix when Lagrange multipliers method is used.
1176 d_residue_d_dofc(disp_range, disp_range) +=
1177 (trans_LR * transposed(trans_LR)) * penalty_factor;
1178 }
1179 if (!constraint_has_penalty) {
1180 trans_LR *= lagrange_mult_scale;
1181 b2linalg::Matrix<T, b2linalg::Mcompressed_col> LR = transposed(trans_LR);
1182 d_residue_d_dofc(lagrange_mult_range, disp_range) += LR;
1183 if (!MATRIX_FORMAT::symmetric) {
1184 d_residue_d_dofc(disp_range, lagrange_mult_range) += trans_LR;
1185 }
1186 }
1187 }
1188 if (!d_residue_d_time.is_null()) {
1189 d_residue_d_time(disp_range) = trans_d_r_d_dof * d_fe_d_time;
1190 }
1191 }
1192
1193 const b2linalg::Index& get_dof_field() {
1194 if (dof_field.empty()
1195 && mrbc->get_size(false, b2linalg::Vector<T, b2linalg::Vdense>::null, 1) != 0) {
1196 b2linalg::Matrix<T, b2linalg::Mcompressed_col> trans_d_r_d_dof;
1197 mrbc->get_nonlinear_value(
1198 b2linalg::Vector<double, b2linalg::Vdense_constref>::null, 0, false,
1199 b2linalg::Vector<T, b2linalg::Vdense_ref>::null, trans_d_r_d_dof,
1200 b2linalg::Matrix<T, b2linalg::Mrectangle_ref>::null, 0);
1201 dof_field.resize(trans_d_r_d_dof.size1());
1202 const b2linalg::Index& df = model->get_domain().get_dof_field().second;
1203 size_t* colind;
1204 size_t* rowind;
1205 T* value;
1206 const size_t s = trans_d_r_d_dof.size2();
1207 trans_d_r_d_dof.get_values(colind, rowind, value);
1208 const size_t* rowind_ptr = rowind;
1209 for (size_t j = 0; j != s; ++j) {
1210 const size_t* rowind_end = rowind + colind[j + 1];
1211 for (; rowind_ptr != rowind_end; ++rowind_ptr) {
1212 dof_field[*rowind_ptr] = std::max(dof_field[*rowind_ptr], df[j]);
1213 }
1214 }
1215 }
1216 return dof_field;
1217 }
1218
1219 typedef ObjectTypeComplete<
1220 DynamicOrderNResidueFunction, typename DynamicResidueFunction<T, MATRIX_FORMAT>::type_t>
1221 type_t;
1222 static type_t type;
1223
1224protected:
1225 logging::Logger& logger;
1226 size_t number_of_dof;
1227 size_t number_of_dof_red;
1228 Model* model;
1229 Domain* domain;
1230 mrbc_t* mrbc;
1231 ebc_t* ebc;
1232 nbc_t* nbc;
1233 SolutionDynamicNonlinear<T>* solution;
1234 b2linalg::Vector<T, b2linalg::Vdense> fi;
1235 b2linalg::Vector<T, b2linalg::Vdense> nbc_sol;
1236 b2linalg::Index ebc_to_remove;
1237 SolverUtils<T, MATRIX_FORMAT> solver_utile;
1238 bool constraint_has_penalty;
1239 double penalty_factor;
1240 double lagrange_mult_scale;
1241 bool rcfo_only_spc;
1242 b2linalg::Index dof_field;
1243};
1244
1245template <typename T, typename MATRIX_FORMAT>
1246typename DynamicOrderNResidueFunction<T, MATRIX_FORMAT>::type_t
1247 DynamicOrderNResidueFunction<T, MATRIX_FORMAT>::type(
1248 "DynamicOrderNResidueFunction", type_name<T, MATRIX_FORMAT>(),
1249 StringList(
1250 "ORDER_ONE_DYNAMIC_RESIDUE_FUNCTION", "ORDER_TWO_DYNAMIC_RESIDUE_FUNCTION",
1251 "ORDER_N_DYNAMIC_RESIDUE_FUNCTION"),
1252 b2000::solver::module, &DynamicResidueFunction<T, MATRIX_FORMAT>::type);
1253
1254template <typename T, typename MATRIX_FORMAT = b2linalg::Msymmetric>
1255class DynamicOrderTwoRayleighDampingResidueFunction
1256 : public DynamicOrderNResidueFunction<T, MATRIX_FORMAT> {
1257public:
1258 typedef DynamicOrderNResidueFunction<T, MATRIX_FORMAT> parent;
1259
1260 void init(Model& model, const AttributeList& attribute) {
1261 parent::init(model, attribute);
1262 if (parent::get_order() != 2) {
1263 Exception() << "Rayleigh damping can only be used for models "
1264 "of order 2"
1265 << THROW;
1266 }
1267 solver_utile.init(
1268 rayleigh_damping_alpha, rayleigh_damping_beta, *parent::model, parent::logger);
1269 }
1270
1271 void set_case(Case* case_) {
1272 rayleigh_damping_alpha = case_->get_double("RAYLEIGH_ALPHA");
1273 rayleigh_damping_beta = case_->get_double("RAYLEIGH_BETA");
1274 }
1275
1276 void save_solution(
1277 b2linalg::Matrix<T, b2linalg::Mrectangle>& dof_red, const double time, const double dtime,
1278 const RTable& attributes, bool end_of_stage = false, bool unconverged_subcycle = false) {
1279 const size_t mrbc_size =
1280 parent::mrbc->get_size(false, b2linalg::Vector<T, b2linalg::Vdense>::null, 1);
1281 b2linalg::Interval disp_range(0, mrbc_size);
1282 b2linalg::Matrix<T, b2linalg::Mrectangle> dof(parent::number_of_dof, dof_red.size2());
1283 b2linalg::Matrix<T, b2linalg::Mcompressed_col> trans_d_r_d_dof;
1284 b2linalg::Matrix<T, b2linalg::Mrectangle> d_r_d_time(
1285 parent::number_of_dof, dof_red.size2() - 1);
1286 parent::mrbc->get_nonlinear_value(
1287 dof_red[0](disp_range), time, true, dof[0], trans_d_r_d_dof, d_r_d_time, 0);
1288 for (size_t j = 1; j < dof.size2(); ++j) {
1289 dof[j] = transposed(trans_d_r_d_dof) * dof_red[j](disp_range);
1290 dof[j] += d_r_d_time[j - 1];
1291 }
1292 if (parent::rcfo_only_spc) {
1293 // recompute fi using only the spc contribution
1294
1295 b2linalg::Matrix<T, b2linalg::Mcompressed_col> trans_d_l_d_dof;
1296 b2linalg::Vector<T> l(parent::ebc->get_size(false));
1297 parent::ebc->get_nonlinear_value(
1298 dof[0], time, false, l, trans_d_l_d_dof,
1299 b2linalg::Vector<T, b2linalg::Vdense_ref>::null, 0);
1300 if (!parent::ebc_to_remove.empty()) {
1301 l.remove(parent::ebc_to_remove);
1302 trans_d_l_d_dof.remove_col(parent::ebc_to_remove);
1303 }
1304
1305 if (parent::constraint_has_penalty) {
1306 l *= -parent::penalty_factor;
1307 } else {
1308 b2linalg::Interval lagrange_mult_range(mrbc_size, dof_red.size1());
1309 if (MATRIX_FORMAT::symmetric && parent::penalty_factor != 0) {
1310 // We apply penalty method with a small factor to obtain a definite positive
1311 // matrix.
1312 l *= -parent::penalty_factor;
1313 l += -parent::lagrange_mult_scale * dof_red[0](lagrange_mult_range);
1314 } else {
1315 l = -parent::lagrange_mult_scale * dof_red[0](lagrange_mult_range);
1316 }
1317 }
1318 for (size_t j = 0; j != trans_d_l_d_dof.size2(); ++j) {
1319 if (trans_d_l_d_dof.get_nb_nonzero(j) != 1) { l[j] = 0; }
1320 }
1321 parent::fi = trans_d_l_d_dof * l;
1322 }
1323 if (unconverged_subcycle) {
1324 parent::solution->set_solution(
1325 dof, time, parent::nbc_sol, parent::fi, unconverged_subcycle);
1326 return;
1327 }
1328 parent::solution->new_cycle(end_of_stage);
1329 parent::solution->set_solution(dof, time, parent::nbc_sol, parent::fi);
1330
1331 typename SolutionDynamicNonlinear<T>::gradient_container gc =
1332 parent::solution->get_gradient_container();
1333 parent::domain->set_dof(dof);
1334 EquilibriumSolution es(1, 0);
1335 double energy[3];
1336 solver_utile.get_elements_values(
1337 dof, time, dtime, es, false, trans_d_r_d_dof, d_r_d_time[0],
1338 b2linalg::Vector<double, b2linalg::Vdense>::null, *parent::model,
1339 b2linalg::Vector<T, b2linalg::Vdense>::null,
1340 b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>::null,
1341 b2linalg::Vector<T, b2linalg::Vdense>::null, gc.get(), 0, parent::logger, 0, energy);
1342 parent::solution->set_value("RAYLEIGH_KINEMATIC_ENERGY", energy[0]);
1343 parent::solution->set_value("RAYLEIGH_RATE_DAMPING_DISSIPATED_ENERGY", energy[1]);
1344 parent::solution->set_value("RAYLEIGH_DAMPING_DISSIPATED_ENERGY", energy[2]);
1345 std::string domain_state_id = parent::solution->get_domain_state_id();
1346 if (!domain_state_id.empty()) { parent::domain->save_state(domain_state_id); }
1347 parent::solution->terminate_cycle(time, dtime, attributes);
1348 if (end_of_stage) { parent::solution->terminate_stage(true); }
1349 }
1350
1351 void get_residue(
1352 const b2linalg::Matrix<T, b2linalg::Mrectangle>& dof_red, const double time,
1353 const double delta_time, const EquilibriumSolution equilibrium_solution,
1354 b2linalg::Vector<T, b2linalg::Vdense>& residue,
1355 const b2linalg::Vector<double, b2linalg::Vdense>& d_residue_d_dof_coeff,
1356 b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>& d_residue_d_dofc,
1357 b2linalg::Vector<T, b2linalg::Vdense>& d_residue_d_time, SolverHints* solver_hints,
1358 CorrectionTerminationTest<T>* termination_test = 0) {
1359 // logging::get_logger("solver.static_nonlinear.residue_function")
1360 b2linalg::Interval disp_range(
1361 0, parent::mrbc->get_size(false, b2linalg::Vector<T, b2linalg::Vdense>::null, 1));
1362 b2linalg::Interval lagrange_mult_range(
1363 parent::mrbc->get_size(false, b2linalg::Vector<T, b2linalg::Vdense>::null, 1),
1364 parent::number_of_dof_red);
1365
1366 if (!residue.is_null()) { residue.set_zero(); }
1367 if (!d_residue_d_dofc.is_null()) { d_residue_d_dofc.set_zero_same_struct(); }
1368
1369 b2linalg::Matrix<T, b2linalg::Mcompressed_col> trans_d_r_d_dof;
1370 b2linalg::Matrix<T, b2linalg::Mrectangle> d_r_d_time(
1371 parent::domain->get_number_of_dof(), dof_red.size2() - 1);
1372 b2linalg::Matrix<T, b2linalg::Mrectangle> dof(parent::number_of_dof, dof_red.size2());
1373 parent::mrbc->get_nonlinear_value(
1374 dof_red[0](disp_range), time, equilibrium_solution, dof[0], trans_d_r_d_dof,
1375 d_r_d_time, solver_hints);
1376 for (size_t j = 1; j < dof.size2(); ++j) {
1377 dof[j] = transposed(trans_d_r_d_dof) * dof_red[j](disp_range);
1378 dof[j] += d_r_d_time[j - 1];
1379 }
1380 b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible> d_fe_d_dof;
1381 b2linalg::Vector<T, b2linalg::Vdense> d_fe_d_time(parent::domain->get_number_of_dof());
1382 parent::nbc_sol.set_zero();
1383 parent::nbc->get_nonlinear_value(
1384 dof[0], time, equilibrium_solution,
1385 residue.is_null() ? b2linalg::Vector<T, b2linalg::Vdense_ref>::null
1386 : b2linalg::Vector<T, b2linalg::Vdense_ref>(parent::nbc_sol),
1387 d_residue_d_dofc.is_null()
1388 ? b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>::null
1389 : d_fe_d_dof,
1390 d_residue_d_time.is_null() ? b2linalg::Vector<T, b2linalg::Vdense_ref>::null
1391 : b2linalg::Vector<T, b2linalg::Vdense_ref>(d_fe_d_time),
1392 solver_hints);
1393 if (termination_test && !residue.is_null()) { termination_test->add_flux(parent::nbc_sol); }
1394 parent::fi = -parent::nbc_sol;
1395 d_fe_d_time = -d_fe_d_time;
1396 if (d_fe_d_dof.size1() != 0) {
1397 d_residue_d_dofc += -(trans_d_r_d_dof * d_fe_d_dof * transposed(trans_d_r_d_dof));
1398 }
1399
1400 b2linalg::Matrix<T, b2linalg::Mcompressed_col> trans_d_l_d_dof;
1401 b2linalg::Vector<T> l;
1402
1403 if (parent::constraint_has_penalty) {
1404 if (!residue.is_null()) { l.resize(parent::ebc->get_size(false)); }
1405 parent::ebc->get_nonlinear_value(
1406 dof[0], time, equilibrium_solution,
1407 (residue.is_null() ? b2linalg::Vector<T, b2linalg::Vdense_ref>::null
1408 : b2linalg::Vector<T, b2linalg::Vdense_ref>(l)),
1409 trans_d_l_d_dof, b2linalg::Vector<T, b2linalg::Vdense_ref>::null, solver_hints);
1410 } else {
1411 if (parent::ebc_to_remove.empty()) {
1412 parent::ebc->get_nonlinear_value(
1413 dof[0], time, equilibrium_solution,
1414 residue.is_null() ? b2linalg::Vector<T, b2linalg::Vdense_ref>::null
1415 : residue(lagrange_mult_range),
1416 trans_d_l_d_dof,
1417 d_residue_d_time.is_null() ? b2linalg::Vector<T, b2linalg::Vdense_ref>::null
1418 : d_residue_d_time(lagrange_mult_range),
1419 solver_hints);
1420 } else {
1421 b2linalg::Vector<T> l(residue.is_null() ? 0 : parent::ebc->get_size(false));
1422 parent::ebc->get_nonlinear_value(
1423 dof[0], time, equilibrium_solution,
1424 (residue.is_null() ? b2linalg::Vector<T, b2linalg::Vdense_ref>::null
1425 : b2linalg::Vector<T, b2linalg::Vdense_ref>(l)),
1426 trans_d_l_d_dof,
1427 (d_residue_d_time.is_null() ? b2linalg::Vector<T, b2linalg::Vdense_ref>::null
1428 : d_residue_d_time(lagrange_mult_range)),
1429 solver_hints);
1430 if (!residue.is_null()) {
1431 l.remove(parent::ebc_to_remove);
1432 residue(lagrange_mult_range) = l;
1433 }
1434 trans_d_l_d_dof.remove_col(parent::ebc_to_remove);
1435 }
1436 if (!d_residue_d_time.is_null()) {
1437 d_residue_d_time(lagrange_mult_range) +=
1438 transposed(trans_d_l_d_dof) * d_r_d_time[0];
1439 }
1440 }
1441
1442 if (!d_residue_d_time.is_null() && d_fe_d_dof.size1()) {
1443 d_fe_d_time -= d_fe_d_dof * d_r_d_time[0];
1444 }
1445
1446 solver_utile.get_elements_values(
1447 dof, time, delta_time, equilibrium_solution, false, trans_d_r_d_dof, d_r_d_time[0],
1448 d_residue_d_dof_coeff, *parent::model, parent::fi, d_residue_d_dofc,
1449 d_residue_d_time.is_null() ? b2linalg::Vector<T, b2linalg::Vdense_ref>::null
1450 : b2linalg::Vector<T, b2linalg::Vdense_ref>(d_fe_d_time),
1451 0, solver_hints, parent::logger, termination_test);
1452
1453 if (!residue.is_null()) {
1454 b2linalg::Vector<T> tmp = parent::fi;
1455 if (parent::constraint_has_penalty) {
1456 tmp += (trans_d_l_d_dof * l) * parent::penalty_factor;
1457 } else {
1458 if (MATRIX_FORMAT::symmetric && parent::penalty_factor != 0) {
1459 // We apply penalty method with a small factor to obtain a definite positive
1460 // matrix.
1461 b2linalg::Vector<T> tmp1;
1462 tmp1 = residue(lagrange_mult_range) * parent::penalty_factor;
1463 tmp1 += dof_red[0](lagrange_mult_range) * parent::lagrange_mult_scale;
1464 tmp += trans_d_l_d_dof * tmp1;
1465 } else {
1466 tmp += parent::lagrange_mult_scale
1467 * (trans_d_l_d_dof * dof_red[0](lagrange_mult_range));
1468 }
1469 residue(lagrange_mult_range) *= parent::lagrange_mult_scale;
1470 }
1471 residue(disp_range) = trans_d_r_d_dof * tmp;
1472 }
1473 if (!d_residue_d_dofc.is_null()) {
1474 b2linalg::Matrix<T, b2linalg::Mcompressed_col> trans_LR;
1475 trans_LR = trans_d_r_d_dof * trans_d_l_d_dof;
1476 if (parent::penalty_factor != 0
1477 && (parent::constraint_has_penalty || MATRIX_FORMAT::symmetric)) {
1478 // We apply penalty method but with a small factor to
1479 // obtain a definite positive matrix when Lagrange multipliers method is used.
1480 d_residue_d_dofc(disp_range, disp_range) +=
1481 (trans_LR * transposed(trans_LR)) * parent::penalty_factor;
1482 }
1483 if (!parent::constraint_has_penalty) {
1484 trans_LR *= parent::lagrange_mult_scale;
1485 b2linalg::Matrix<T, b2linalg::Mcompressed_col> LR = transposed(trans_LR);
1486 d_residue_d_dofc(lagrange_mult_range, disp_range) += LR;
1487 if (!MATRIX_FORMAT::symmetric) {
1488 d_residue_d_dofc(disp_range, lagrange_mult_range) += trans_LR;
1489 }
1490 }
1491 }
1492 if (!d_residue_d_time.is_null()) {
1493 d_residue_d_time(disp_range) = trans_d_r_d_dof * d_fe_d_time;
1494 }
1495 }
1496
1497 typedef ObjectTypeComplete<
1498 DynamicOrderTwoRayleighDampingResidueFunction,
1499 typename DynamicResidueFunction<T, MATRIX_FORMAT>::type_t>
1500 type_t;
1501 static type_t type;
1502
1503private:
1504 double rayleigh_damping_alpha;
1505 double rayleigh_damping_beta;
1506 SolverUtilsOrderTwoRayleighDamping<T, MATRIX_FORMAT> solver_utile;
1507};
1508
1509template <typename T, typename MATRIX_FORMAT>
1510typename DynamicOrderTwoRayleighDampingResidueFunction<T, MATRIX_FORMAT>::type_t
1511 DynamicOrderTwoRayleighDampingResidueFunction<T, MATRIX_FORMAT>::type(
1512 "DynamicOrderTwoRayleighDampingResidueFunction", type_name<T, MATRIX_FORMAT>(),
1513 StringList("RAYLEIGH_DAMPING_DYNAMIC_RESIDUE_FUNCTION"), b2000::solver::module,
1514 &DynamicResidueFunction<T, MATRIX_FORMAT>::type);
1515
1516template <typename T, typename MATRIX_FORMAT = b2linalg::Msymmetric>
1517class QuasistaticViscousDampingResidueFunction
1518 : public DynamicOrderNResidueFunction<T, MATRIX_FORMAT> {
1519public:
1520 typedef DynamicOrderNResidueFunction<T, MATRIX_FORMAT> parent;
1521
1522 void init(Model& model, const AttributeList& attribute) {
1523 parent::init(model, attribute);
1524 if (parent::get_order() != 2) {
1525 Exception() << "Quasi-static viscous damping can only be "
1526 "used for models of order 2"
1527 << THROW;
1528 }
1529 solver_utile.init(
1530 rayleigh_damping_alpha, rayleigh_damping_beta, *parent::model, parent::logger);
1531 }
1532
1533 void set_case(Case* case_) {
1534 const char* alpha_s = "RAYLEIGH_ALPHA";
1535 const char* beta_s = "RAYLEIGH_BETA";
1536 if (case_->has_key(alpha_s) && case_->has_key(beta_s)) {
1537 rayleigh_damping_alpha = case_->get_double("RAYLEIGH_ALPHA");
1538 rayleigh_damping_beta = case_->get_double("RAYLEIGH_BETA");
1539 dissipated_energy_fraction = 0;
1540 } else {
1541 dissipated_energy_fraction = case_->get_double("DISSIPATED_ENERGY_FRACTION", 1.e-4);
1542 rayleigh_damping_alpha = 1e-10; // we set a very small dissipation for the first step
1543 rayleigh_damping_beta = 0;
1544 }
1545 }
1546
1547 int get_order() const { return 1; }
1548
1549 void save_solution(
1550 b2linalg::Matrix<T, b2linalg::Mrectangle>& dof_red, const double time, const double dtime,
1551 const RTable& attributes, bool end_of_stage = false, bool unconverged_subcycle = false) {
1552 const size_t mrbc_size =
1553 parent::mrbc->get_size(false, b2linalg::Vector<T, b2linalg::Vdense>::null, 1);
1554 b2linalg::Interval disp_range(0, mrbc_size);
1555 b2linalg::Matrix<T, b2linalg::Mrectangle> dof(parent::number_of_dof, dof_red.size2());
1556 b2linalg::Matrix<T, b2linalg::Mcompressed_col> trans_d_r_d_dof;
1557 b2linalg::Matrix<T, b2linalg::Mrectangle> d_r_d_time(
1558 parent::number_of_dof, dof_red.size2() - 1);
1559 parent::mrbc->get_nonlinear_value(
1560 dof_red[0](disp_range), time, true, dof[0], trans_d_r_d_dof, d_r_d_time, 0);
1561 for (size_t j = 1; j < dof.size2(); ++j) {
1562 dof[j] = transposed(trans_d_r_d_dof) * dof_red[j](disp_range);
1563 dof[j] += d_r_d_time[j - 1];
1564 }
1565 if (parent::rcfo_only_spc) {
1566 // recompute fi using only the spc contribution
1567
1568 b2linalg::Matrix<T, b2linalg::Mcompressed_col> trans_d_l_d_dof;
1569 b2linalg::Vector<T> l(parent::ebc->get_size(false));
1570 parent::ebc->get_nonlinear_value(
1571 dof[0], time, false, l, trans_d_l_d_dof,
1572 b2linalg::Vector<T, b2linalg::Vdense_ref>::null, 0);
1573 if (!parent::ebc_to_remove.empty()) {
1574 l.remove(parent::ebc_to_remove);
1575 trans_d_l_d_dof.remove_col(parent::ebc_to_remove);
1576 }
1577
1578 if (parent::constraint_has_penalty) {
1579 l *= -parent::penalty_factor;
1580 } else {
1581 b2linalg::Interval lagrange_mult_range(mrbc_size, dof_red.size1());
1582 if (MATRIX_FORMAT::symmetric && parent::penalty_factor != 0) {
1583 // We apply penalty method with a small factor to obtain a definite positive
1584 // matrix.
1585 l *= -parent::penalty_factor;
1586 l += -parent::lagrange_mult_scale * dof_red[0](lagrange_mult_range);
1587 } else {
1588 l = -parent::lagrange_mult_scale * dof_red[0](lagrange_mult_range);
1589 }
1590 }
1591 for (size_t j = 0; j != trans_d_l_d_dof.size2(); ++j) {
1592 if (trans_d_l_d_dof.get_nb_nonzero(j) != 1) { l[j] = 0; }
1593 }
1594 parent::fi = trans_d_l_d_dof * l;
1595 }
1596 if (unconverged_subcycle) {
1597 parent::solution->set_solution(
1598 dof, time, parent::nbc_sol, parent::fi, unconverged_subcycle);
1599 return;
1600 }
1601 parent::solution->new_cycle(end_of_stage);
1602 parent::solution->set_solution(dof, time, parent::nbc_sol, parent::fi);
1603
1604 typename SolutionDynamicNonlinear<T>::gradient_container gc =
1605 parent::solution->get_gradient_container();
1606 parent::domain->set_dof(dof);
1607 EquilibriumSolution es(1, 0);
1608 double energy[3];
1609 solver_utile.get_elements_values(
1610 dof, time, dtime, es, false, trans_d_r_d_dof, d_r_d_time[0],
1611 b2linalg::Vector<double, b2linalg::Vdense>::null, *parent::model,
1612 b2linalg::Vector<T, b2linalg::Vdense>::null,
1613 b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>::null,
1614 b2linalg::Vector<T, b2linalg::Vdense>::null, gc.get(), 0, parent::logger, 0, energy);
1615 parent::solution->set_value("RAYLEIGH_KINEMATIC_ENERGY", energy[0]);
1616 parent::solution->set_value("RAYLEIGH_RATE_DAMPING_DISSIPATED_ENERGY", energy[1]);
1617 parent::solution->set_value("RAYLEIGH_DAMPING_DISSIPATED_ENERGY", energy[2]);
1618 if (dissipated_energy_fraction != 0) {
1619 double res = 0;
1620 for (size_t i = 0; i != dof.size1(); ++i) {
1621 // res += std::max(norm(dof(i, 0) * parent::nbc_sol[i]), norm(dof(i, 0) *
1622 // parent::fi[i]));
1623 res += dof(i, 0) * (parent::nbc_sol[i] + parent::fi[i]);
1624 }
1625 double sfactor = dissipated_energy_fraction * res / (energy[1] * dtime);
1626 parent::logger << logging::debug << "scale factor for the stabilised method "
1627 << sfactor * 1e-10 << logging::LOGFLUSH;
1628 solver_utile.scale_raylay_damping(sfactor);
1629 dissipated_energy_fraction = 0;
1630 }
1631 std::string domain_state_id = parent::solution->get_domain_state_id();
1632 if (!domain_state_id.empty()) { parent::domain->save_state(domain_state_id); }
1633 parent::solution->terminate_cycle(time, dtime, attributes);
1634 if (end_of_stage) { parent::solution->terminate_stage(true); }
1635 }
1636
1637 void get_residue(
1638 const b2linalg::Matrix<T, b2linalg::Mrectangle>& dof_red, const double time,
1639 const double delta_time, const EquilibriumSolution equilibrium_solution,
1640 b2linalg::Vector<T, b2linalg::Vdense>& residue,
1641 const b2linalg::Vector<double, b2linalg::Vdense>& d_residue_d_dof_coeff,
1642 b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>& d_residue_d_dofc,
1643 b2linalg::Vector<T, b2linalg::Vdense>& d_residue_d_time, SolverHints* solver_hints,
1644 CorrectionTerminationTest<T>* termination_test = 0) {
1645 // logging::get_logger("solver.static_nonlinear.residue_function")
1646 b2linalg::Interval disp_range(
1647 0, parent::mrbc->get_size(false, b2linalg::Vector<T, b2linalg::Vdense>::null, 1));
1648 b2linalg::Interval lagrange_mult_range(
1649 parent::mrbc->get_size(false, b2linalg::Vector<T, b2linalg::Vdense>::null, 1),
1650 parent::number_of_dof_red);
1651
1652 if (!residue.is_null()) { residue.set_zero(); }
1653 if (!d_residue_d_dofc.is_null()) { d_residue_d_dofc.set_zero_same_struct(); }
1654
1655 b2linalg::Matrix<T, b2linalg::Mcompressed_col> trans_d_r_d_dof;
1656 b2linalg::Matrix<T, b2linalg::Mrectangle> d_r_d_time(
1657 parent::domain->get_number_of_dof(), dof_red.size2() - 1);
1658 b2linalg::Matrix<T, b2linalg::Mrectangle> dof(parent::number_of_dof, dof_red.size2());
1659 parent::mrbc->get_nonlinear_value(
1660 dof_red[0](disp_range), time, equilibrium_solution, dof[0], trans_d_r_d_dof,
1661 d_r_d_time, solver_hints);
1662 for (size_t j = 1; j < dof.size2(); ++j) {
1663 dof[j] = transposed(trans_d_r_d_dof) * dof_red[j](disp_range);
1664 dof[j] += d_r_d_time[j - 1];
1665 }
1666 b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible> d_fe_d_dof;
1667 b2linalg::Vector<T, b2linalg::Vdense> d_fe_d_time(parent::domain->get_number_of_dof());
1668 parent::nbc_sol.set_zero();
1669 parent::nbc->get_nonlinear_value(
1670 dof[0], time, equilibrium_solution,
1671 residue.is_null() ? b2linalg::Vector<T, b2linalg::Vdense_ref>::null
1672 : b2linalg::Vector<T, b2linalg::Vdense_ref>(parent::nbc_sol),
1673 d_residue_d_dofc.is_null()
1674 ? b2linalg::Matrix<T, typename MATRIX_FORMAT::sparse_inversible>::null
1675 : d_fe_d_dof,
1676 d_residue_d_time.is_null() ? b2linalg::Vector<T, b2linalg::Vdense_ref>::null
1677 : b2linalg::Vector<T, b2linalg::Vdense_ref>(d_fe_d_time),
1678 solver_hints);
1679 if (termination_test && !residue.is_null()) { termination_test->add_flux(parent::nbc_sol); }
1680 parent::fi = -parent::nbc_sol;
1681 d_fe_d_time = -d_fe_d_time;
1682 if (d_fe_d_dof.size1() != 0) {
1683 d_residue_d_dofc += -(trans_d_r_d_dof * d_fe_d_dof * transposed(trans_d_r_d_dof));
1684 }
1685
1686 b2linalg::Matrix<T, b2linalg::Mcompressed_col> trans_d_l_d_dof;
1687 b2linalg::Vector<T> l;
1688
1689 if (parent::constraint_has_penalty) {
1690 if (!residue.is_null()) { l.resize(parent::ebc->get_size(false)); }
1691 parent::ebc->get_nonlinear_value(
1692 dof[0], time, equilibrium_solution,
1693 (residue.is_null() ? b2linalg::Vector<T, b2linalg::Vdense_ref>::null
1694 : b2linalg::Vector<T, b2linalg::Vdense_ref>(l)),
1695 trans_d_l_d_dof, b2linalg::Vector<T, b2linalg::Vdense_ref>::null, solver_hints);
1696 } else {
1697 if (parent::ebc_to_remove.empty()) {
1698 parent::ebc->get_nonlinear_value(
1699 dof[0], time, equilibrium_solution,
1700 residue.is_null() ? b2linalg::Vector<T, b2linalg::Vdense_ref>::null
1701 : residue(lagrange_mult_range),
1702 trans_d_l_d_dof,
1703 d_residue_d_time.is_null() ? b2linalg::Vector<T, b2linalg::Vdense_ref>::null
1704 : d_residue_d_time(lagrange_mult_range),
1705 solver_hints);
1706 } else {
1707 b2linalg::Vector<T> l(residue.is_null() ? 0 : parent::ebc->get_size(false));
1708 parent::ebc->get_nonlinear_value(
1709 dof[0], time, equilibrium_solution,
1710 (residue.is_null() ? b2linalg::Vector<T, b2linalg::Vdense_ref>::null
1711 : b2linalg::Vector<T, b2linalg::Vdense_ref>(l)),
1712 trans_d_l_d_dof,
1713 (d_residue_d_time.is_null() ? b2linalg::Vector<T, b2linalg::Vdense_ref>::null
1714 : d_residue_d_time(lagrange_mult_range)),
1715 solver_hints);
1716 if (!residue.is_null()) {
1717 l.remove(parent::ebc_to_remove);
1718 residue(lagrange_mult_range) = l;
1719 }
1720 trans_d_l_d_dof.remove_col(parent::ebc_to_remove);
1721 }
1722 if (!d_residue_d_time.is_null()) {
1723 d_residue_d_time(lagrange_mult_range) +=
1724 transposed(trans_d_l_d_dof) * d_r_d_time[0];
1725 }
1726 }
1727
1728 if (!d_residue_d_time.is_null() && d_fe_d_dof.size1()) {
1729 d_fe_d_time -= d_fe_d_dof * d_r_d_time[0];
1730 }
1731
1732 solver_utile.get_elements_values(
1733 dof, time, delta_time, equilibrium_solution, false, trans_d_r_d_dof, d_r_d_time[0],
1734 d_residue_d_dof_coeff, *parent::model, parent::fi, d_residue_d_dofc,
1735 d_residue_d_time.is_null() ? b2linalg::Vector<T, b2linalg::Vdense_ref>::null
1736 : b2linalg::Vector<T, b2linalg::Vdense_ref>(d_fe_d_time),
1737 0, solver_hints, parent::logger, termination_test);
1738
1739 if (!residue.is_null()) {
1740 b2linalg::Vector<T> tmp = parent::fi;
1741 if (parent::constraint_has_penalty) {
1742 tmp += (trans_d_l_d_dof * l) * parent::penalty_factor;
1743 } else {
1744 if (MATRIX_FORMAT::symmetric && parent::penalty_factor != 0) {
1745 // We apply penalty method with a small factor to obtain a definite positive
1746 // matrix.
1747 b2linalg::Vector<T> tmp1;
1748 tmp1 = residue(lagrange_mult_range) * parent::penalty_factor;
1749 tmp1 += dof_red[0](lagrange_mult_range) * parent::lagrange_mult_scale;
1750 tmp += trans_d_l_d_dof * tmp1;
1751 } else {
1752 tmp += parent::lagrange_mult_scale
1753 * (trans_d_l_d_dof * dof_red[0](lagrange_mult_range));
1754 }
1755 residue(lagrange_mult_range) *= parent::lagrange_mult_scale;
1756 }
1757 residue(disp_range) = trans_d_r_d_dof * tmp;
1758 }
1759 if (!d_residue_d_dofc.is_null()) {
1760 b2linalg::Matrix<T, b2linalg::Mcompressed_col> trans_LR;
1761 trans_LR = trans_d_r_d_dof * trans_d_l_d_dof;
1762 if (parent::penalty_factor != 0
1763 && (parent::constraint_has_penalty || MATRIX_FORMAT::symmetric)) {
1764 // We apply penalty method but with a small factor to
1765 // obtain a definite positive matrix when Lagrange multipliers method is used.
1766 d_residue_d_dofc(disp_range, disp_range) +=
1767 (trans_LR * transposed(trans_LR)) * parent::penalty_factor;
1768 }
1769 if (!parent::constraint_has_penalty) {
1770 trans_LR *= parent::lagrange_mult_scale;
1771 b2linalg::Matrix<T, b2linalg::Mcompressed_col> LR = transposed(trans_LR);
1772 d_residue_d_dofc(lagrange_mult_range, disp_range) += LR;
1773 if (!MATRIX_FORMAT::symmetric) {
1774 d_residue_d_dofc(disp_range, lagrange_mult_range) += trans_LR;
1775 }
1776 }
1777 }
1778 if (!d_residue_d_time.is_null()) {
1779 d_residue_d_time(disp_range) = trans_d_r_d_dof * d_fe_d_time;
1780 }
1781 }
1782
1783 typedef ObjectTypeComplete<
1784 QuasistaticViscousDampingResidueFunction,
1785 typename DynamicResidueFunction<T, MATRIX_FORMAT>::type_t>
1786 type_t;
1787 static type_t type;
1788
1789private:
1790 double rayleigh_damping_alpha;
1791 double rayleigh_damping_beta;
1792 double dissipated_energy_fraction;
1793 SolverUtilsOrderTwoRayleighDamping<T, MATRIX_FORMAT, false> solver_utile;
1794};
1795
1796template <typename T, typename MATRIX_FORMAT>
1797typename QuasistaticViscousDampingResidueFunction<T, MATRIX_FORMAT>::type_t
1798 QuasistaticViscousDampingResidueFunction<T, MATRIX_FORMAT>::type(
1799 "QuasistaticViscousDampingResidueFunction", type_name<T, MATRIX_FORMAT>(),
1800 StringList(
1801 "QUASISTATIC_VISCOUS_DAMPING_DYNAMIC_RESIDUE_FUNCTION",
1802 "ARTIFICIAL_DAMPING_DYNAMIC_RESIDUE_FUNCTION"),
1803 b2000::solver::module, &DynamicResidueFunction<T, MATRIX_FORMAT>::type);
1804
1805}} // namespace b2000::solver
1806
1807#endif // B2DYNAMIC_NONLINEAR_UTILE_H_
b2boundary_condition.H
b2000::details::pow
csda< T > pow(const csda< T > &a, const csda< T > &b)
Power from two csda numbers.
Definition b2csda.H:378
b2000::details::norm
csda< T > norm(const csda< T > &a)
Definition b2csda.H:343
THROW
#define THROW
Definition b2exception.H:198
b2solver.H
Interface to C++ representations of FE solvers.
b2000::SolverHints::diverge
@ diverge
Definition b2solver.H:183
b2000::SolverHints::terminate_analysis
@ terminate_analysis
Definition b2solver.H:194
b2000::logging::get_logger
Logger & get_logger(const std::string &logger_name="")
Definition b2logging.H:829
b2000
Contains the base classes for implementing Finite Elements.
Definition b2boundary_condition.H:32