-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathtest_userfunctions.cpp
More file actions
114 lines (105 loc) · 6.52 KB
/
Copy pathtest_userfunctions.cpp
File metadata and controls
114 lines (105 loc) · 6.52 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
#include <cstdio>
#include <cmath>
#include <limits>
#include "UserFunctions.H"
int main() {
int fail = 0;
// ---- eigensystem: SWE quasilinear A (n=2), u=2, c=sqrt(g h)=3
// A = [[0,1],[c^2-u^2, 2u]] -> eigenvalues u-c=-1, u+c=5
const double u = 2.0, c = 3.0;
double a[4] = {0.0, 1.0, c*c - u*u, 2*u};
double lam[2], R[4], L[4];
for (int i = 0; i < 2; ++i) lam[i] = eigensystem(i, a[0], a[1], a[2], a[3]);
for (int i = 0; i < 4; ++i) R[i] = eigensystem(2 + i, a[0], a[1], a[2], a[3]);
for (int i = 0; i < 4; ++i) L[i] = eigensystem(6 + i, a[0], a[1], a[2], a[3]);
printf("eig lambda = %.6f %.6f (expect -1, 5 in some order)\n", lam[0], lam[1]);
double lo = std::min(lam[0], lam[1]), hi = std::max(lam[0], lam[1]);
if (std::fabs(lo - (u - c)) > 1e-10 || std::fabs(hi - (u + c)) > 1e-10) { printf(" FAIL eigenvalues\n"); fail++; }
// reconstruct A = R * diag(lam) * L -> must recover the input exactly
double err = 0.0;
for (int i = 0; i < 2; ++i) for (int j = 0; j < 2; ++j) {
double s = 0.0;
for (int k = 0; k < 2; ++k) s += R[i*2+k] * lam[k] * L[k*2+j];
err = std::max(err, std::fabs(s - a[i*2+j]));
}
printf("||R diag(lam) L - A||_inf = %.3e (consistency of the SAME eigenbasis)\n", err);
if (err > 1e-10) { printf(" FAIL reconstruction\n"); fail++; }
// L must be the inverse of R
double ierr = 0.0;
for (int i = 0; i < 2; ++i) for (int j = 0; j < 2; ++j) {
double s = 0.0;
for (int k = 0; k < 2; ++k) s += L[i*2+k] * R[k*2+j];
ierr = std::max(ierr, std::fabs(s - (i == j ? 1.0 : 0.0)));
}
printf("||L R - I||_inf = %.3e\n", ierr);
if (ierr > 1e-10) { printf(" FAIL L != R^-1\n"); fail++; }
// Roe dissipation |A| = R |Lam| L must be well defined & symmetric-consistent
double absA[4];
for (int i = 0; i < 2; ++i) for (int j = 0; j < 2; ++j) {
double s = 0.0;
for (int k = 0; k < 2; ++k) s += R[i*2+k] * std::fabs(lam[k]) * L[k*2+j];
absA[i*2+j] = s;
}
printf("|A| = [[%.4f, %.4f], [%.4f, %.4f]]\n", absA[0], absA[1], absA[2], absA[3]);
// ---- eigenvalues: the lambda-only kernel must agree with eigensystem's lambda
double ev[2];
for (int i = 0; i < 2; ++i) ev[i] = eigenvalues(i, a[0], a[1], a[2], a[3]);
printf("eigenvalues = %.6f %.6f (expect -1, 5)\n", ev[0], ev[1]);
double evlo = std::min(ev[0], ev[1]), evhi = std::max(ev[0], ev[1]);
if (std::fabs(evlo - (u - c)) > 1e-10 || std::fabs(evhi - (u + c)) > 1e-10) { printf(" FAIL eigenvalues\n"); fail++; }
// must match eigensystem's spectrum as a SET (order is not guaranteed)
if (std::fabs(evlo - lo) > 1e-10 || std::fabs(evhi - hi) > 1e-10) { printf(" FAIL eigenvalues != eigensystem spectrum\n"); fail++; }
// max|lambda| is what local_max_abs_eigenvalue actually consumes
double mx = std::max(std::fabs(ev[0]), std::fabs(ev[1]));
printf("max|lambda| = %.6f (expect 5 = |u|+c)\n", mx);
if (std::fabs(mx - 5.0) > 1e-10) { printf(" FAIL max|lambda|\n"); fail++; }
// cache invalidation on a different matrix
double e2 = eigenvalues(0, 2.0, 0.0, 0.0, 3.0);
double e3 = eigenvalues(1, 2.0, 0.0, 0.0, 3.0);
printf("eigenvalues(diag(2,3)) = %.6f %.6f (expect 2, 3 -> cache invalidates)\n", e2, e3);
if (std::fabs(std::min(e2,e3) - 2.0) > 1e-10 || std::fabs(std::max(e2,e3) - 3.0) > 1e-10) { printf(" FAIL eigenvalues cache\n"); fail++; }
// ---- solve: A x = b, n=3, known answer
// A = [[2,1,0],[1,3,1],[0,1,4]], b = [3,9,13] -> x = [1,1,3]
double x0 = solve(0, 2.,1.,0., 1.,3.,1., 0.,1.,4., 3.,7.,13.);
double x1 = solve(1, 2.,1.,0., 1.,3.,1., 0.,1.,4., 3.,7.,13.);
double x2 = solve(2, 2.,1.,0., 1.,3.,1., 0.,1.,4., 3.,7.,13.);
printf("solve x = %.10f %.10f %.10f (expect 1, 1, 3)\n", x0, x1, x2);
if (std::fabs(x0-1) > 1e-10 || std::fabs(x1-1) > 1e-10 || std::fabs(x2-3) > 1e-10) { printf(" FAIL solve\n"); fail++; }
// ---- cache invalidation: a DIFFERENT matrix must not return the stale answer
double y0 = solve(0, 1.,0.,0., 0.,1.,0., 0.,0.,1., 7.,8.,9.);
double y1 = solve(1, 1.,0.,0., 0.,1.,0., 0.,0.,1., 7.,8.,9.);
printf("solve(I, [7,8,9]) = %.6f %.6f (expect 7, 8 -> cache invalidates)\n", y0, y1);
if (std::fabs(y0-7) > 1e-10 || std::fabs(y1-8) > 1e-10) { printf(" FAIL cache invalidation\n"); fail++; }
// ---- REQ-168 ADDENDUM (1): inf-tolerance. MOOD feeds non-finite candidate
// states by design; the kernel must report +inf, not garbage. (Eigen does
// NOT raise -- it returns NoConvergence and leaves the values unset.)
{
const double INF = std::numeric_limits<double>::infinity();
const double QNAN = std::nan("");
double ei = eigenvalues(0, 0.0, 1.0, INF, 4.0);
double en = eigenvalues(0, 0.0, 1.0, QNAN, 4.0);
printf("eigenvalues(inf matrix) = %f eigenvalues(nan matrix) = %f (expect +inf both)\n", ei, en);
if (!std::isinf(ei) || ei < 0) { printf(" FAIL inf-tolerance (inf)\n"); fail++; }
if (!std::isinf(en) || en < 0) { printf(" FAIL inf-tolerance (nan)\n"); fail++; }
// eigensystem: lambda=+inf, R=L=I -> |A| = R|Lam|L = inf*I, never NaN
double lam0 = eigensystem(0, 0.0, 1.0, INF, 4.0);
double R00 = eigensystem(2, 0.0, 1.0, INF, 4.0);
double R01 = eigensystem(3, 0.0, 1.0, INF, 4.0);
double L00 = eigensystem(6, 0.0, 1.0, INF, 4.0);
printf("eigensystem(inf): lam0=%f R=[%g %g ..] L00=%g (expect inf, I, I)\n", lam0, R00, R01, L00);
if (!std::isinf(lam0)) { printf(" FAIL eigensystem inf lambda\n"); fail++; }
if (R00 != 1.0 || R01 != 0.0 || L00 != 1.0) { printf(" FAIL eigensystem inf R/L != I\n"); fail++; }
// and the |A| product must be inf, NOT NaN (0*inf would poison the row)
double absA00 = R00 * std::fabs(lam0) * L00;
printf(" |A|[0,0] = %f isnan=%d (must be inf, never NaN)\n", absA00, (int)std::isnan(absA00));
if (std::isnan(absA00)) { printf(" FAIL |A| is NaN\n"); fail++; }
// a finite matrix must still work AFTER a non-finite one (cache key)
double back = eigenvalues(0, 2.0, 0.0, 0.0, 3.0);
double back1 = eigenvalues(1, 2.0, 0.0, 0.0, 3.0);
printf(" recovery after non-finite: %.4f %.4f (expect 2, 3)\n", back, back1);
if (std::fabs(std::min(back,back1)-2.0) > 1e-10 || std::fabs(std::max(back,back1)-3.0) > 1e-10) {
printf(" FAIL recovery after non-finite\n"); fail++; }
}
printf(fail ? "\nRESULT: %d FAILURES\n" : "\nRESULT: all pass\n", fail);
return fail;
}