ryujin 2.1.1 revision 5b63c078534504e9d75978bd539bf7ffb02e1157
riemann_solver.template.h
Go to the documentation of this file.
1//
2// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
3// Copyright (C) 2020 - 2024 by the ryujin authors
4//
5
6#pragma once
7
8#include <compile_time_options.h>
9
10#include "riemann_solver.h"
11
12#include <newton.h>
13#include <simd.h>
14
15// #define DEBUG_RIEMANN_SOLVER
16
17namespace ryujin
18{
19 namespace Euler
20 {
21 template <int dim, typename Number>
22 DEAL_II_ALWAYS_INLINE inline Number
24 const Number p_star) const
25 {
26 const auto view = hyperbolic_system.view<dim, Number>();
27 const auto &gamma = view.gamma();
28
29 const auto &[rho, u, p, a] = riemann_data;
30
31 const Number Az = ScalarNumber(2.) / (rho * (gamma + Number(1.)));
32 const Number Bz =
33 (gamma - ScalarNumber(1.)) / (gamma + ScalarNumber(1.)) * p;
34 const Number radicand = Az / (p_star + Bz);
35 const Number true_value = (p_star - p) * std::sqrt(radicand);
36
37 const auto exponent =
38 ScalarNumber(0.5) * (gamma - ScalarNumber(1.)) / gamma;
39 const Number factor = ryujin::pow(p_star / p, exponent) - Number(1.);
40 const auto false_value =
41 ScalarNumber(2.) * a * factor / (gamma - ScalarNumber(1.));
42
43 return dealii::compare_and_apply_mask<
44 dealii::SIMDComparison::greater_than_or_equal>(
45 p_star, p, true_value, false_value);
46 }
47
48
49 template <int dim, typename Number>
50 DEAL_II_ALWAYS_INLINE inline Number
52 const Number &p_star) const
53 {
54 const auto view = hyperbolic_system.view<dim, Number>();
55
57 const auto &gamma = view.gamma();
58 const auto &gamma_inverse = view.gamma_inverse();
59 const auto &gamma_minus_one_inverse = view.gamma_minus_one_inverse();
60 const auto &gamma_plus_one_inverse = view.gamma_plus_one_inverse();
61
62 const auto &[rho, u, p, a] = riemann_data;
63
64 const Number radicand_inverse = ScalarNumber(0.5) * rho *
65 ((gamma + ScalarNumber(1.)) * p_star +
66 (gamma - ScalarNumber(1.)) * p);
67 const Number denominator =
68 (p_star + (gamma - ScalarNumber(1.)) * gamma_plus_one_inverse * p);
69 const Number true_value =
70 (denominator - ScalarNumber(0.5) * (p_star - p)) /
71 (denominator * std::sqrt(radicand_inverse));
72
73 const auto exponent =
74 (ScalarNumber(-1.) - gamma) * ScalarNumber(0.5) * gamma_inverse;
75 const Number factor = (gamma - ScalarNumber(1.)) * ScalarNumber(0.5) *
76 gamma_inverse * ryujin::pow(p_star / p, exponent) /
77 p;
78 const auto false_value =
79 factor * ScalarNumber(2.) * a * gamma_minus_one_inverse;
80
81 return dealii::compare_and_apply_mask<
82 dealii::SIMDComparison::greater_than_or_equal>(
83 p_star, p, true_value, false_value);
84 }
85
86
87 template <int dim, typename Number>
88 DEAL_II_ALWAYS_INLINE inline Number
90 const primitive_type &riemann_data_j,
91 const Number p_in) const
92 {
93 const Number &u_i = riemann_data_i[1];
94 const Number &u_j = riemann_data_j[1];
95
96 return f(riemann_data_i, p_in) + f(riemann_data_j, p_in) + u_j - u_i;
97 }
98
99
100 template <int dim, typename Number>
101 DEAL_II_ALWAYS_INLINE inline Number
103 const primitive_type &riemann_data_j,
104 const Number &p) const
105 {
106 return df(riemann_data_i, p) + df(riemann_data_j, p);
107 }
108
109
110 /*
111 * The approximate Riemann solver is based on a function phi(p) that is
112 * montone increasing in p, concave down and whose (weak) third
113 * derivative is non-negative and locally bounded [1, p. 912]. Because we
114 * actually do not perform any iteration for computing our wavespeed
115 * estimate we can get away by only implementing a specialized variant of
116 * the phi function that computes phi(p_max). It inlines the
117 * implementation of the "f" function and eliminates all unnecessary
118 * branches in "f".
119 *
120 * Cost: 0x pow, 2x division, 2x sqrt
121 */
122 template <int dim, typename Number>
123 DEAL_II_ALWAYS_INLINE inline Number
125 const primitive_type &riemann_data_i,
126 const primitive_type &riemann_data_j) const
127 {
128 const auto view = hyperbolic_system.view<dim, Number>();
129 const auto &gamma = view.gamma();
130
131 const auto &[rho_i, u_i, p_i, a_i] = riemann_data_i;
132 const auto &[rho_j, u_j, p_j, a_j] = riemann_data_j;
133
134 const Number p_max = std::max(p_i, p_j);
135
136 const Number radicand_inverse_i = ScalarNumber(0.5) * rho_i *
137 ((gamma + ScalarNumber(1.)) * p_max +
138 (gamma - ScalarNumber(1.)) * p_i);
139
140 const Number value_i = (p_max - p_i) / std::sqrt(radicand_inverse_i);
141
142 const Number radicand_inverse_j = ScalarNumber(0.5) * rho_j *
143 ((gamma + ScalarNumber(1.)) * p_max +
144 (gamma - ScalarNumber(1.)) * p_j);
145
146 const Number value_j = (p_max - p_j) / std::sqrt(radicand_inverse_j);
147
148 return value_i + value_j + u_j - u_i;
149 }
150
151
152 /*
153 * Next we construct approximations for the two extreme wave speeds of
154 * the Riemann fan [1, p. 912, (3.7) + (3.8)] and compute an upper bound
155 * lambda_max of the maximal wave speed:
156 */
157
158
159 /*
160 * see [1], page 912, (3.7)
161 *
162 * Cost: 0x pow, 1x division, 1x sqrt
163 */
164 template <int dim, typename Number>
165 DEAL_II_ALWAYS_INLINE inline Number
167 const primitive_type &riemann_data, const Number p_star) const
168 {
169 const auto view = hyperbolic_system.view<dim, Number>();
170 const auto &gamma = view.gamma();
171 const auto &gamma_inverse = view.gamma_inverse();
172 const auto factor =
173 (gamma + ScalarNumber(1.0)) * ScalarNumber(0.5) * gamma_inverse;
174
175 const auto &[rho, u, p, a] = riemann_data;
176 const auto inv_p = ScalarNumber(1.0) / p;
177
178 const Number tmp = positive_part((p_star - p) * inv_p);
179
180 return u - a * std::sqrt(ScalarNumber(1.0) + factor * tmp);
181 }
182
183
184 /*
185 * see [1], page 912, (3.8)
186 *
187 * Cost: 0x pow, 1x division, 1x sqrt
188 */
189 template <int dim, typename Number>
190 DEAL_II_ALWAYS_INLINE inline Number
192 const primitive_type &primitive_state, const Number p_star) const
193 {
194 const auto view = hyperbolic_system.view<dim, Number>();
195 const auto &gamma = view.gamma();
196 const auto &gamma_inverse = view.gamma_inverse();
197 const Number factor =
198 (gamma + ScalarNumber(1.0)) * ScalarNumber(0.5) * gamma_inverse;
199
200 const auto &[rho, u, p, a] = primitive_state;
201 const auto inv_p = ScalarNumber(1.0) / p;
202
203 const Number tmp = positive_part((p_star - p) * inv_p);
204 return u + a * std::sqrt(Number(1.0) + factor * tmp);
205 }
206
207
217 template <int dim, typename Number>
218 DEAL_II_ALWAYS_INLINE inline std::array<Number, 2>
220 const std::array<Number, 4> &riemann_data_i,
221 const std::array<Number, 4> &riemann_data_j,
222 const Number p_1,
223 const Number p_2) const
224 {
225 const Number nu_11 = lambda1_minus(riemann_data_i, p_2 /*SIC!*/);
226 const Number nu_12 = lambda1_minus(riemann_data_i, p_1 /*SIC!*/);
227
228 const Number nu_31 = lambda3_plus(riemann_data_j, p_1);
229 const Number nu_32 = lambda3_plus(riemann_data_j, p_2);
230
231 const Number lambda_max =
232 std::max(positive_part(nu_32), negative_part(nu_11));
233
234 const Number gap =
235 std::max(std::abs(nu_32 - nu_31), std::abs(nu_12 - nu_11));
236
237 return {{gap, lambda_max}};
238 }
239
240
241 /*
242 * For two given primitive states <code>riemann_data_i</code> and
243 * <code>riemann_data_j</code>, and a guess p_2, compute an upper bound
244 * for lambda.
245 *
246 * This is the same lambda_max as computed by compute_gap. The function
247 * simply avoids a number of unnecessary computations (in case we do
248 * not need to know the gap).
249 *
250 * Cost: 0x pow, 2x division, 2x sqrt
251 */
252 template <int dim, typename Number>
253 DEAL_II_ALWAYS_INLINE inline Number
255 const primitive_type &riemann_data_i,
256 const primitive_type &riemann_data_j,
257 const Number p_star) const
258 {
259 const Number nu_11 = lambda1_minus(riemann_data_i, p_star);
260 const Number nu_32 = lambda3_plus(riemann_data_j, p_star);
261
262 return std::max(positive_part(nu_32), negative_part(nu_11));
263 }
264
265
266 /*
267 * Two-rarefaction approximation to p_star computed for two primitive
268 * states <code>riemann_data_i</code> and <code>riemann_data_j</code>.
269 *
270 * See [1], page 914, (4.3)
271 *
272 * Cost: 2x pow, 2x division, 0x sqrt
273 */
274 template <int dim, typename Number>
275 DEAL_II_ALWAYS_INLINE inline Number
277 const primitive_type &riemann_data_i,
278 const primitive_type &riemann_data_j) const
279 {
280 const auto view = hyperbolic_system.view<dim, Number>();
281 const auto &gamma = view.gamma();
282 const auto &gamma_inverse = view.gamma_inverse();
283 const auto &gamma_minus_one_inverse = view.gamma_minus_one_inverse();
284
285 const auto &[rho_i, u_i, p_i, a_i] = riemann_data_i;
286 const auto &[rho_j, u_j, p_j, a_j] = riemann_data_j;
287 const auto inv_p_j = ScalarNumber(1.) / p_j;
288
289 /*
290 * Nota bene (cf. [1, (4.3)]):
291 * a_Z^0 * sqrt(1 - b * rho_Z) = a_Z * (1 - b * rho_Z)
292 * We have computed a_Z already, so we are simply going to use this
293 * identity below:
294 */
295
296 const auto factor = (gamma - ScalarNumber(1.)) * ScalarNumber(0.5);
297
298 /*
299 * Nota bene (cf. [1, (3.6)]: The condition "numerator > 0" is the
300 * well-known non-vacuum condition. In case we encounter numerator <= 0
301 * then p_star = 0 is the correct pressure to compute the wave speed.
302 * Therefore, all we have to do is to take the positive part of the
303 * expression:
304 */
305
306 const Number numerator = positive_part(a_i + a_j - factor * (u_j - u_i));
307 const Number denominator =
308 a_i * ryujin::pow(p_i * inv_p_j, -factor * gamma_inverse) + a_j;
309
310 const auto exponent = ScalarNumber(2.0) * gamma * gamma_minus_one_inverse;
311
312 const auto p_1_tilde =
313 p_j * ryujin::pow(numerator / denominator, exponent);
314
315#ifdef DEBUG_RIEMANN_SOLVER
316 std::cout << "p_star_two_rarefaction = " << p_1_tilde << std::endl;
317#endif
318 return p_1_tilde;
319 }
320
321
322 /*
323 * Failsafe approximation to p_star computed for two primitive
324 * states <code>riemann_data_i</code> and <code>riemann_data_j</code>.
325 *
326 * See [1], page 914, (4.3)
327 *
328 * Cost: 2x pow, 2x division, 0x sqrt
329 */
330 template <int dim, typename Number>
331 DEAL_II_ALWAYS_INLINE inline Number
333 const primitive_type &riemann_data_i,
334 const primitive_type &riemann_data_j) const
335 {
336 const auto view = hyperbolic_system.view<dim, Number>();
337 const auto &gamma = view.gamma();
338
339 const auto &[rho_i, u_i, p_i, a_i] = riemann_data_i;
340 const auto &[rho_j, u_j, p_j, a_j] = riemann_data_j;
341
342 /*
343 * Compute (5.11) formula for \tilde p_2^\ast:
344 *
345 * Cost: 0x pow, 3x division, 3x sqrt
346 */
347
348 const Number p_max = std::max(p_i, p_j);
349
350 Number radicand_i = ScalarNumber(2.) * p_max;
351 radicand_i /=
352 rho_i * ((gamma + Number(1.)) * p_max + (gamma - Number(1.)) * p_i);
353
354 const Number x_i = std::sqrt(radicand_i);
355
356 Number radicand_j = ScalarNumber(2.) * p_max;
357 radicand_j /=
358 rho_j * ((gamma + Number(1.)) * p_max + (gamma - Number(1.)) * p_j);
359
360 const Number x_j = std::sqrt(radicand_j);
361
362 const Number a = x_i + x_j;
363 const Number b = u_j - u_i;
364 const Number c = -p_i * x_i - p_j * x_j;
365
366 const Number base = (-b + std::sqrt(b * b - ScalarNumber(4.) * a * c)) /
367 (ScalarNumber(2.) * a);
368 const Number p_2_tilde = base * base;
369
370#ifdef DEBUG_RIEMANN_SOLVER
371 std::cout << "p_star_failsafe = " << p_2_tilde << std::endl;
372#endif
373 return p_2_tilde;
374 }
375
376
377 template <int dim, typename Number>
378 DEAL_II_ALWAYS_INLINE inline auto
380 const state_type &U, const dealii::Tensor<1, dim, Number> &n_ij) const
382 {
383 const auto view = hyperbolic_system.view<dim, Number>();
384
385 const auto rho = view.density(U);
386 const auto rho_inverse = Number(1.0) / rho;
387
388 const auto m = view.momentum(U);
389 const auto proj_m = n_ij * m;
390 const auto perp = m - proj_m * n_ij;
391
392 const auto E =
393 view.total_energy(U) - Number(0.5) * perp.norm_square() * rho_inverse;
394
395 using state_type_1d =
397 const auto view_1d = hyperbolic_system.view<1, Number>();
398
399 const auto state = state_type_1d{{rho, proj_m, E}};
400 const auto p = view_1d.pressure(state);
401 const auto a = view_1d.speed_of_sound(state);
402 return {{rho, proj_m * rho_inverse, p, a}};
403 }
404
405
406 template <int dim, typename Number>
408 const primitive_type &riemann_data_i,
409 const primitive_type &riemann_data_j) const
410 {
411 /*
412 * For exactly solving the Riemann problem we need to start with a
413 * good upper and lower bound, p_1 <= p_star <= p_2, for finding
414 * phi(p_star) == 0. This implies that we have to ensure that
415 * phi(p_2) >= 0 and phi(p_1) <= 0.
416 *
417 * Instead of solving the Riemann problem exactly, however we will
418 * simply use the upper bound p_2 (with p_2 >= p_star) to compute
419 * lambda_max and return the estimate.
420 *
421 * We will use three candidates, p_min, p_max and the two rarefaction
422 * approximation p_star_tilde. We have (up to round-off errors) that
423 * phi(p_star_tilde) >= 0. So this is a safe upper bound, it might
424 * just be too large.
425 *
426 * Depending on the sign of phi(p_max) we select the following ranges:
427 *
428 * phi(p_max) < 0:
429 * p_1 <- p_max and p_2 <- p_star_tilde
430 *
431 * phi(p_max) >= 0:
432 * p_1 <- p_min and p_2 <- min(p_max, p_star_tilde)
433 *
434 * Nota bene:
435 *
436 * - The special case phi(p_max) == 0 as discussed in [1] is already
437 * contained in the second condition.
438 *
439 * - In principle, we would have to treat the case phi(p_min) > 0 as
440 * well. This corresponds to two expansion waves and a good
441 * estimate for the wavespeed is obtained by simply computing
442 * lambda_max with p_2 = 0.
443 *
444 * However, it turns out that numerically in this case we will
445 * have
446 *
447 * 0 < p_star <= p_star_tilde <= p_min <= p_max.
448 *
449 * So it is sufficient to end up with p_2 = p_star_tilde (!!) to
450 * compute the exact same wave speed as for p_2 = 0.
451 *
452 * Note: If for some reason p_star should be computed exactly,
453 * then p_1 has to be set to zero. This can be done efficiently by
454 * simply checking for p_2 < p_1 and setting p_1 <- 0 if
455 * necessary.
456 */
457
458 const auto &[rho_i, u_i, p_i, a_i] = riemann_data_i;
459 const auto &[rho_j, u_j, p_j, a_j] = riemann_data_j;
460
461#ifdef DEBUG_RIEMANN_SOLVER
462 std::cout << "rho_left: " << rho_i << std::endl;
463 std::cout << "u_left: " << u_i << std::endl;
464 std::cout << "p_left: " << p_i << std::endl;
465 std::cout << "a_left: " << a_i << std::endl;
466 std::cout << "rho_right: " << rho_j << std::endl;
467 std::cout << "u_right: " << u_j << std::endl;
468 std::cout << "p_right: " << p_j << std::endl;
469 std::cout << "a_right: " << a_j << std::endl;
470#endif
471
472 const Number p_max = std::max(p_i, p_j);
473
474 const Number rarefaction =
475 p_star_two_rarefaction(riemann_data_i, riemann_data_j);
476 const Number failsafe = p_star_failsafe(riemann_data_i, riemann_data_j);
477 const Number p_star_tilde = std::min(rarefaction, failsafe);
478
479 const Number phi_p_max = phi_of_p_max(riemann_data_i, riemann_data_j);
480
481 Number p_2 =
482 dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
483 phi_p_max,
484 Number(0.),
485 p_star_tilde,
486 std::min(p_max, p_star_tilde));
487
488#ifdef DEBUG_RIEMANN_SOLVER
489 std::cout << " p^*_tilde = " << p_2 << "\n";
490 std::cout << " phi(p_*_t) = "
491 << phi(riemann_data_i, riemann_data_j, p_2) << std::endl;
492#endif
493
494 /*
495 * If we do no Newton iteration, cut it short:
496 */
497
498 if (parameters.newton_max_iterations() == 0) {
499 const auto lambda_max =
500 compute_lambda(riemann_data_i, riemann_data_j, p_2);
501
502#ifdef DEBUG_RIEMANN_SOLVER
503 std::cout << "-> lambda_max = " << lambda_max << std::endl;
504#endif
505 return lambda_max;
506 }
507
508 /*
509 * Compute p_1 and ensure that p_1 < p_2. If we hit a case with two
510 * expansions we might indeed have that p_star_tilde < p_1. Set p_1 =
511 * p_2 in this case.
512 */
513
514 const Number p_min = std::min(riemann_data_i[2], riemann_data_j[2]);
515
516 Number p_1 =
517 dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
518 phi_p_max, Number(0.), p_max, p_min);
519
520 p_1 = dealii::compare_and_apply_mask<
521 dealii::SIMDComparison::less_than_or_equal>(p_1, p_2, p_1, p_2);
522
523 /*
524 * Step 2: Perform quadratic Newton iteration.
525 *
526 * See [1], p. 915f (4.8) and (4.9)
527 */
528
529 auto [gap, lambda_max] =
530 compute_gap(riemann_data_i, riemann_data_j, p_1, p_2);
531
532#ifdef DEBUG_RIEMANN_SOLVER
533 std::cout << std::fixed << std::setprecision(16);
534 std::cout << "p_1: (start) " << p_1 << std::endl;
535 std::cout << "p_2: (start) " << p_2 << std::endl;
536 std::cout << "gap: (start) " << gap << std::endl;
537 std::cout << "l_m: (start) " << lambda_max << std::endl;
538#endif
539
540 for (unsigned int i = 0; i < parameters.newton_max_iterations(); ++i) {
541
542 /* We accept our current guess if we reach the tolerance... */
543 const Number tolerance(parameters.newton_tolerance());
544 if (std::max(Number(0.), gap - tolerance) == Number(0.)) {
545#ifdef DEBUG_RIEMANN_SOLVER
546 std::cout << "converged after " << i << " iterations." << std::endl;
547#endif
548 break;
549 }
550
551 // FIXME: Fuse these computations:
552 const Number phi_p_1 = phi(riemann_data_i, riemann_data_j, p_1);
553 const Number phi_p_2 = phi(riemann_data_i, riemann_data_j, p_2);
554 const Number dphi_p_1 = dphi(riemann_data_i, riemann_data_j, p_1);
555 const Number dphi_p_2 = dphi(riemann_data_i, riemann_data_j, p_2);
556
557 quadratic_newton_step(p_1, p_2, phi_p_1, phi_p_2, dphi_p_1, dphi_p_2);
558
559 /* Update lambda_max and gap: */
560 auto [gap_new, lambda_max_new] =
561 compute_gap(riemann_data_i, riemann_data_j, p_1, p_2);
562 gap = gap_new;
563 lambda_max = lambda_max_new;
564
565#ifdef DEBUG_RIEMANN_SOLVER
566 std::cout << "phi_p_1: " << phi_p_1 << std::endl;
567 std::cout << "phi_p_2: " << phi_p_2 << std::endl;
568 std::cout << "dphi_p_1: " << dphi_p_1 << std::endl;
569 std::cout << "dphi_p_2: " << dphi_p_2 << std::endl;
570 std::cout << "p_1: ( " << i << " ) " << p_1 << std::endl;
571 std::cout << "p_2: ( " << i << " ) " << p_2 << std::endl;
572 std::cout << "gap: " << gap << std::endl;
573 std::cout << "l_m: " << lambda_max << std::endl;
574#endif
575 }
576
577#ifdef DEBUG_RIEMANN_SOLVER
578 std::cout << "-> lambda_max = " << lambda_max << std::endl;
579#endif
580
581 return lambda_max;
582 }
583
584
585 template <int dim, typename Number>
586 DEAL_II_ALWAYS_INLINE inline Number RiemannSolver<dim, Number>::compute(
587 const state_type &U_i,
588 const state_type &U_j,
589 const unsigned int /*i*/,
590 const unsigned int * /*js*/,
591 const dealii::Tensor<1, dim, Number> &n_ij) const
592 {
593 const auto riemann_data_i = riemann_data_from_state(U_i, n_ij);
594 const auto riemann_data_j = riemann_data_from_state(U_j, n_ij);
595
596 return compute(riemann_data_i, riemann_data_j);
597 }
598
599 } // namespace Euler
600} // namespace ryujin
dealii::Tensor< 1, problem_dimension, Number > state_type
Number phi_of_p_max(const primitive_type &riemann_data_i, const primitive_type &riemann_data_j) const
Number compute_lambda(const primitive_type &riemann_data_i, const primitive_type &riemann_data_j, const Number p_star) const
typename View::ScalarNumber ScalarNumber
Number df(const primitive_type &riemann_data, const Number &p_star) const
Number lambda1_minus(const primitive_type &riemann_data, const Number p_star) const
Number phi(const primitive_type &riemann_data_i, const primitive_type &riemann_data_j, const Number p_in) const
Number p_star_two_rarefaction(const primitive_type &riemann_data_i, const primitive_type &riemann_data_j) const
primitive_type riemann_data_from_state(const state_type &U, const dealii::Tensor< 1, dim, Number > &n_ij) const
std::array< Number, 2 > compute_gap(const primitive_type &riemann_data_i, const primitive_type &riemann_data_j, const Number p_1, const Number p_2) const
typename View::state_type state_type
Number compute(const primitive_type &riemann_data_i, const primitive_type &riemann_data_j) const
std::array< Number, riemann_data_size > primitive_type
Number p_star_failsafe(const primitive_type &riemann_data_i, const primitive_type &riemann_data_j) const
Number dphi(const primitive_type &riemann_data_i, const primitive_type &riemann_data_j, const Number &p) const
Number f(const primitive_type &riemann_data, const Number p_star) const
Number lambda3_plus(const primitive_type &primitive_state, const Number p_star) const
T pow(const T x, const T b)
DEAL_II_ALWAYS_INLINE Number negative_part(const Number number)
Definition: simd.h:124
DEAL_II_ALWAYS_INLINE Number positive_part(const Number number)
Definition: simd.h:112
DEAL_II_ALWAYS_INLINE void quadratic_newton_step(Number &p_1, Number &p_2, const Number phi_p_1, const Number phi_p_2, const Number dphi_p_1, const Number dphi_p_2, const Number sign=Number(1.0))
Definition: newton.h:39