v0.16.0
Loading...
Searching...
No Matches
matrix_function.cpp
Go to the documentation of this file.
1/**
2 * @file matrix_function.cpp
3 * @example mofem/atom_tests/matrix_function.cpp
4 * @brief Test and example for matrix function
5 *
6 * For reference see \cite miehe2001algorithms
7 *
8 */
9
10//#define FTENSOR_DEBUG
11#include <algorithm>
12
13#include <FTensor.hpp>
14
15#include <MoFEM.hpp>
16using namespace MoFEM;
17
18#include <MatrixFunction.hpp>
19
20FTensor::Index<'i', 3> i;
21FTensor::Index<'j', 3> j;
22FTensor::Index<'k', 3> k;
23FTensor::Index<'l', 3> l;
24
25template <typename T1, typename T2, int DIM>
26void diff_ddg(T1 &t_1, T2 &t_2, const FTensor::Number<DIM> &) {
27 constexpr double eps = 1e-4;
28 for (int ii = 0; ii != DIM; ++ii)
29 for (int jj = 0; jj != DIM; ++jj)
30 for (int kk = 0; kk != DIM; ++kk)
31 for (int ll = 0; ll != DIM; ++ll) {
32
33 if (std::abs(t_1(ii, jj, kk, ll) - t_2(ii, jj, kk, ll)) > eps)
34 MOFEM_LOG("ATOM_TEST", Sev::error)
35 << "Error " << ii << " " << jj << " " << kk << " " << ll << " "
36 << t_1(ii, jj, kk, ll) << " " << t_2(ii, jj, kk, ll);
37 }
38
39 for (int ii = 0; ii != DIM; ++ii)
40 for (int jj = 0; jj != DIM; ++jj)
41 for (int kk = 0; kk != DIM; ++kk)
42 for (int ll = 0; ll != DIM; ++ll)
43 t_1(ii, jj, kk, ll) -= t_2(ii, jj, kk, ll);
44};
45
46template <typename T1, typename T2, int DIM>
47auto get_diff_matrix2(T1 &t_a, T2 &t_d, const FTensor::Number<DIM> &) {
48 constexpr auto t_kd = FTensor::Kronecker_Delta<double>();
49 FTensor::Tensor4<double, DIM, DIM, DIM, DIM> t_d_a;
50
51 FTensor::Index<'i', DIM> i;
52 FTensor::Index<'j', DIM> j;
53 FTensor::Index<'k', DIM> k;
54 FTensor::Index<'l', DIM> l;
55
56 t_d_a(i, j, k, l) = 0;
57
58 for (int ii = 0; ii != DIM; ++ii)
59 for (int jj = 0; jj != DIM; ++jj)
60 for (int kk = 0; kk != DIM; ++kk)
61 for (int ll = 0; ll != DIM; ++ll)
62 for (int zz = 0; zz != DIM; ++zz) {
63
64 auto diff = [&](auto ii, auto jj, auto kk, auto ll, int zz) {
65 return
66
67 t_a(ii, zz) * t_kd(zz, kk) * t_kd(jj, ll)
68
69 +
70
71 t_kd(ii, kk) * t_kd(zz, ll) * t_a(zz, jj);
72 };
73
74 t_d_a(ii, jj, kk, ll) +=
75 (diff(ii, jj, kk, ll, zz) + diff(ii, jj, ll, kk, zz)) / 2.;
76 }
77
78 diff_ddg(t_d_a, t_d, FTensor::Number<DIM>());
79
80 return t_d_a;
81};
82
83template <typename T1, typename T2, int DIM>
84auto get_diff2_matrix2(T1 &t_s, T2 &t_dd, const FTensor::Number<DIM> &) {
85 constexpr auto t_kd = FTensor::Kronecker_Delta<double>();
86 FTensor::Tensor4<double, DIM, DIM, DIM, DIM> t_dd_a;
87
88 FTensor::Index<'i', DIM> i;
89 FTensor::Index<'j', DIM> j;
90 FTensor::Index<'k', DIM> k;
91 FTensor::Index<'l', DIM> l;
92
93 t_dd_a(i, j, k, l) = 0;
94
95 for (int ii = 0; ii != DIM; ++ii)
96 for (int jj = 0; jj != DIM; ++jj)
97 for (int kk = 0; kk != DIM; ++kk)
98 for (int ll = 0; ll != DIM; ++ll)
99 for (int mm = 0; mm != DIM; ++mm)
100 for (int nn = 0; nn != DIM; ++nn)
101 for (int zz = 0; zz != DIM; ++zz) {
102
103 auto diff = [&](auto ii, auto jj, auto kk, auto ll, int mm,
104 int nn, int zz) {
105 return
106
107 t_s(ii, jj) *
108 (t_kd(ii, mm) * t_kd(zz, nn) * t_kd(zz, kk) * t_kd(jj, ll)
109
110 +
111
112 t_kd(ii, kk) * t_kd(zz, ll) * t_kd(zz, mm) *
113 t_kd(jj, nn));
114 };
115
116 t_dd_a(kk, ll, mm, nn) += (
117
118 diff(ii, jj, kk, ll, mm, nn, zz)
119
120 +
121
122 diff(ii, jj, ll, kk, mm, nn, zz)
123
124 +
125
126 diff(ii, jj, kk, ll, nn, mm, zz)
127
128 +
129
130 diff(ii, jj, ll, kk, nn, mm, zz)
131
132 ) /
133 4.;
134 }
135
136 diff_ddg(t_dd_a, t_dd, FTensor::Number<DIM>());
137
138 return t_dd_a;
139};
140
141template <typename T1, int DIM>
142auto get_diff_matrix(T1 &t_d, const FTensor::Number<DIM> &) {
143 constexpr auto t_kd = FTensor::Kronecker_Delta<double>();
144 FTensor::Tensor4<double, DIM, DIM, DIM, DIM> t_d_a;
145
146 FTensor::Index<'i', DIM> i;
147 FTensor::Index<'j', DIM> j;
148 FTensor::Index<'k', DIM> k;
149 FTensor::Index<'l', DIM> l;
150
151 t_d_a(i, j, k, l) = 0;
152
153 for (int ii = 0; ii != DIM; ++ii)
154 for (int jj = 0; jj != DIM; ++jj)
155 for (int kk = 0; kk != DIM; ++kk)
156 for (int ll = 0; ll != DIM; ++ll) {
157
158 auto diff = [&](auto ii, auto jj, auto kk, auto ll) {
159 return t_kd(ii, kk) * t_kd(jj, ll);
160 };
161
162 t_d_a(ii, jj, kk, ll) =
163 (diff(ii, jj, kk, ll) + diff(ii, jj, ll, kk)) / 2.;
164 }
165
166 diff_ddg(t_d_a, t_d, FTensor::Number<3>());
167
168 return t_d_a;
169};
170
171template <typename T, int DIM>
172auto get_norm_t4(T &t, const FTensor::Number<DIM> &) {
173 double r = 0;
174 for (int ii = 0; ii != DIM; ++ii)
175 for (int jj = 0; jj != DIM; ++jj)
176 for (int kk = 0; kk != DIM; ++kk)
177 for (int ll = 0; ll != DIM; ++ll)
178 r += t(ii, jj, kk, ll) * t(ii, jj, kk, ll);
179 return r;
180};
181
182static char help[] = "...\n\n";
183
184int main(int argc, char *argv[]) {
185
186 MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
187
188 auto core_log = logging::core::get();
189 core_log->add_sink(
190 LogManager::createSink(LogManager::getStrmSelf(), "ATOM_TEST"));
191 LogManager::setLog("ATOM_TEST");
192 BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
193 MOFEM_LOG_ATTRIBUTES("ATOM_TEST", LogManager::BitLineID);
194
195 try {
196
197 auto print_ddg = [](auto &t, auto str = "") {
198 constexpr double eps = 1e-6;
199 for (int ii = 0; ii != 3; ++ii)
200 for (int jj = 0; jj != 3; ++jj)
201 for (int kk = 0; kk != 3; ++kk)
202 for (int ll = 0; ll != 3; ++ll) {
203 double v = t(ii, jj, kk, ll);
204 double w = std::abs(v) < eps ? 0 : v;
205 MOFEM_LOG("ATOM_TEST", Sev::noisy)
206 << str << std::fixed << std::setprecision(3) << std::showpos
207 << ii + 1 << " " << jj + 1 << " " << kk + 1 << " " << ll + 1
208 << " : " << w;
209 }
210 };
211
212 auto print_ddg_direction = [](auto &t, auto kk, int ll) {
213 for (int ii = 0; ii != 3; ++ii)
214 for (int jj = 0; jj <= ii; ++jj)
215 MOFEM_LOG("ATOM_TEST", Sev::noisy)
216 << ii + 1 << " " << jj + 1 << " " << kk + 1 << " " << ll + 1
217 << " : " << t(ii, jj, kk, ll);
218 };
219
220 auto print_mat = [](auto &t) {
221 for (int ii = 0; ii != 3; ++ii)
222 for (int jj = 0; jj != 3; ++jj)
223 MOFEM_LOG("ATOM_TEST", Sev::noisy)
224 << ii + 1 << " " << jj + 1 << " : " << t(ii, jj);
225 };
226
227 enum swap { swap12, swap01 };
228 auto run_lapack = [](auto &a, swap s = swap12) {
229 int info;
230 double wkopt;
231 double w[3];
232
233 FTensor::Tensor2<double, 3, 3> t_a{
234
235 a[0], a[1], a[2],
236
237 a[3], a[4], a[5],
238
239 a[6], a[7], a[8]};
240
241 /* Query and allocate the optimal workspace */
242 int lwork = -1;
243 info = lapack_dsyev('V', 'U', 3, a.data(), 3, w, &wkopt, lwork);
244 if (info > 0)
245 THROW_MESSAGE("The algorithm failed to compute eigenvalues.");
246 lwork = (int)wkopt;
247 std::vector<double> work(lwork);
248 /* Solve eigenproblem */
249 info = lapack_dsyev('V', 'U', 3, a.data(), 3, w, &*work.begin(), lwork);
250 if (info > 0)
251 THROW_MESSAGE("The algorithm failed to compute eigenvalues.");
252
253 if (s == swap12) {
254 FTensor::Tensor2<double, 3, 3> t_eig_vec{
255
256 a[0 * 3 + 0], a[0 * 3 + 1], a[0 * 3 + 2],
257
258 a[2 * 3 + 0], a[2 * 3 + 1], a[2 * 3 + 2],
259
260 a[1 * 3 + 0], a[1 * 3 + 1], a[1 * 3 + 2]};
261
262 FTensor::Tensor1<double, 3> t_eig_vals{w[0], w[2], w[1]};
263 return std::make_tuple(t_a, t_eig_vec, t_eig_vals);
264 }
265
266 FTensor::Tensor2<double, 3, 3> t_eig_vec{
267
268 a[1 * 3 + 0], a[1 * 3 + 1], a[1 * 3 + 2],
269
270 a[0 * 3 + 0], a[0 * 3 + 1], a[0 * 3 + 2],
271
272 a[2 * 3 + 0], a[2 * 3 + 1], a[2 * 3 + 2],
273
274 };
275
276 FTensor::Tensor1<double, 3> t_eig_vals{w[1], w[0], w[2]};
277 return std::make_tuple(t_a, t_eig_vec, t_eig_vals);
278 };
279
280 // Test matrix against Mathematica results
281 {
282 FTensor::Tensor2<double, 3, 3> t_A{
283
284 1., 0.1, -0.5,
285
286 0.1, 2., 0.,
287
288 -0.5, 0., 3.};
289
290 FTensor::Tensor2<double, 3, 3> t_N{
291
292 0.969837, -0.0860972, 0.228042,
293
294 0.0790574, 0.996073, 0.0398449,
295
296 -0.230577, -0.0206147, 0.972836};
297
298 FTensor::Tensor1<double, 3> t_L{0.873555, 2.00794, 3.11851};
299
300 FTensor::Tensor2<double, 3, 3> t_S{
301
302 1., 0., 0.,
303
304 0., 1., 0.,
305
306 0., 0., 1.};
307
308 FTensor::Tensor2_symmetric<double, 3> t_S_sym;
309 t_S_sym(i, j) = (t_S(i, j) || t_S(j, i)) / 2.;
310
311 MOFEM_LOG("ATOM_TEST", Sev::verbose) << "Diff A";
312 print_mat(t_A);
313
314 // Test against Mathematica values for the 0th, 1st, and 2nd derivatives
315 {
316 auto f = [](double v) { return exp(v); };
317 auto d_f = [](double v) { return exp(v); };
318 auto dd_f = [](double v) { return exp(v); };
319
320 auto t_b = EigenMatrix::getMat(t_L, t_N, f);
321 MOFEM_LOG("ATOM_TEST", Sev::verbose) << "Reconstruct mat";
322 print_mat(t_b);
323
324 auto t_d = EigenMatrix::getDiffMat(t_L, t_N, f, d_f, 3);
325
326 MOFEM_LOG("ATOM_TEST", Sev::verbose) << "Diff";
327 print_ddg_direction(t_d, 0, 2);
328
329 auto t_dd =
330 EigenMatrix::getDiffDiffMat(t_L, t_N, f, d_f, dd_f, t_S_sym, 3);
331
332 MOFEM_LOG("ATOM_TEST", Sev::verbose) << "Diff Diff";
333 print_ddg_direction(t_dd, 0, 2);
334
335 auto norm2_t_b = t_b(i, j) * t_b(i, j);
336 MOFEM_LOG("ATOM_TEST", Sev::inform) << "norm2_t_b " << norm2_t_b;
337
338 auto norm2_t_d = get_norm_t4(t_d, FTensor::Number<3>());
339 MOFEM_LOG("ATOM_TEST", Sev::inform) << "norm2_t_d " << norm2_t_d;
340
341 auto norm2_t_dd = get_norm_t4(t_dd, FTensor::Number<3>());
342 MOFEM_LOG("ATOM_TEST", Sev::inform) << "norm2_t_dd " << norm2_t_dd;
343
344 constexpr double regression_t_b = 572.543;
345 constexpr double regression_t_d = 859.939;
346 constexpr double regression_t_dd = 859.939;
347
348 constexpr double eps = 1e-2;
349 if (std::abs(norm2_t_b - regression_t_b) > eps)
350 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID, "Wrong t_b");
351 if (std::abs(norm2_t_d - regression_t_d) > eps)
352 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID, "Wrong t_d");
353 if (std::abs(norm2_t_dd - regression_t_dd) > eps)
354 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID, "Wrong t_dd");
355 }
356
357 // Compare with LAPACK-calculated eigenpairs. The `run_lapack()` helper
358 // deliberately reorders eigenpairs for branch testing, so compare spectra
359 // independent of permutation. The matrix function itself should also be
360 // invariant to the sign of the eigenvectors.
361 {
362
363 std::array<double, 9> a{1., 0.1, -0.5,
364
365 0.1, 2., 0.,
366
367 -0.5, 0., 3.};
368
369 auto tuple = run_lapack(a);
370 // auto &t_a = std::get<0>(tuple);
371 auto &t_eig_vec = std::get<1>(tuple);
372 auto &t_eig_vals = std::get<2>(tuple);
373
374 std::array<double, 3> sorted_reference_eigenvalues{
375 t_L(0), t_L(1), t_L(2)};
376 std::array<double, 3> sorted_lapack_eigenvalues{
377 t_eig_vals(0), t_eig_vals(1), t_eig_vals(2)};
378 std::sort(sorted_reference_eigenvalues.begin(),
379 sorted_reference_eigenvalues.end());
380 std::sort(sorted_lapack_eigenvalues.begin(),
381 sorted_lapack_eigenvalues.end());
382
383 double t_eig_val_diff = 0;
384 for (int ii = 0; ii != 3; ++ii)
385 t_eig_val_diff += std::pow(sorted_lapack_eigenvalues[ii] -
386 sorted_reference_eigenvalues[ii],
387 2);
388 MOFEM_LOG("ATOM_TEST", Sev::inform)
389 << "Sorted eigenvalue difference " << t_eig_val_diff;
390
391 auto f = [](double v) { return exp(v); };
392
393 auto t_b = EigenMatrix::getMat(t_L, t_N, f);
394 auto t_c = EigenMatrix::getMat(t_eig_vals, t_eig_vec, f);
395 t_c(i, j) -= t_b(i, j);
396 print_mat(t_c);
397
398 auto norm2_t_c = t_c(i, j) * t_c(i, j);
399 MOFEM_LOG("ATOM_TEST", Sev::inform)
400 << "Reconstruct mat difference with LAPACK eigenvalues and "
401 "vectors "
402 << norm2_t_c;
403
404 constexpr double eps = 1e-8;
405 if (fabs(norm2_t_c) > eps)
406 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
407 "Matrix not reconstructed");
408 }
409
410 constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
411 FTensor::Ddg<double, 3, 3> t_one;
412 t_one(i, j, k, l) = (t_kd(i, k) ^ t_kd(j, l)) / 4.;
413
414 // Test linear function with zero second derivative
415 {
416 auto f = [](double v) { return v; };
417 auto d_f = [](double v) { return 1; };
418 auto dd_f = [](double v) { return 0; };
419
420 constexpr double eps = 1e-10;
421 {
422 auto t_b = EigenMatrix::getMat(t_L, t_N, f);
423 t_b(i, j) -= (t_A(i, j) || t_A(j, i)) / 2;
424 auto norm2_t_b = t_b(i, j) * t_b(i, j);
425 MOFEM_LOG("ATOM_TEST", Sev::inform)
426 << "Result should be matrix itself " << norm2_t_b;
427 if (norm2_t_b > eps)
428 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
429 "This norm should be zero");
430 }
431
432 {
433
434 auto t_d = EigenMatrix::getDiffMat(t_L, t_N, f, d_f, 3);
435 auto t_d_a = get_diff_matrix(t_d, FTensor::Number<3>());
436
437 MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_d_a";
438 print_ddg(t_d_a, "hand ");
439 MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_d";
440 print_ddg(t_d, "code ");
441
442 double nrm2_t_d_a = get_norm_t4(t_d_a, FTensor::Number<3>());
443 MOFEM_LOG("ATOM_TEST", Sev::inform)
444 << "Derivative hand calculation minus code " << nrm2_t_d_a;
445 if (nrm2_t_d_a > eps)
446 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
447 "This norm should be zero");
448 }
449
450 {
451 FTensor::Tensor2<double, 3, 3> t_S{
452
453 1., 0., 0.,
454
455 0., 1., 0.,
456
457 0., 0., 1.};
458
459 FTensor::Tensor2_symmetric<double, 3> t_S_sym;
460 t_S_sym(i, j) = t_S(i, j) || t_S(j, i);
461
462 auto t_dd =
463 EigenMatrix::getDiffDiffMat(t_L, t_N, f, d_f, dd_f, t_S_sym, 3);
464
465 auto norm2_t_dd = get_norm_t4(t_dd, FTensor::Number<3>());
466 MOFEM_LOG("ATOM_TEST", Sev::inform) << "norm2_t_dd " << norm2_t_dd;
467 if (norm2_t_dd > eps)
468 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
469 "This norm should be zero");
470 }
471 }
472
473 // Test quadratic function with constant second derivative
474 {
475 auto f = [](double v) { return v * v; };
476 auto d_f = [](double v) { return 2 * v; };
477 auto dd_f = [](double v) { return 2; };
478
479 constexpr double eps = 1e-9;
480
481 // check if multiplication gives right value
482 {
483 auto t_b = EigenMatrix::getMat(t_L, t_N, f);
484 FTensor::Tensor2<double, 3, 3> t_a;
485 t_a(i, j) = t_b(i, j) - t_A(i, k) * t_A(k, j);
486 print_mat(t_a);
487 auto norm2_t_a = t_a(i, j) * t_a(i, j);
488 MOFEM_LOG("ATOM_TEST", Sev::inform)
489 << "Result should be matrix times matrix " << norm2_t_a;
490 if (norm2_t_a > eps)
491 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
492 "This norm should be zero");
493 }
494
495 // Check first derivative
496 {
497 auto t_d = EigenMatrix::getDiffMat(t_L, t_N, f, d_f, 3);
498 print_ddg_direction(t_d, 0, 2);
499 auto t_d_a = get_diff_matrix2(t_A, t_d, FTensor::Number<3>());
500 MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_d_a";
501 print_ddg(t_d_a, "hand ");
502 MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_d";
503 print_ddg(t_d, "code ");
504 double nrm2_t_d_a = get_norm_t4(t_d_a, FTensor::Number<3>());
505 MOFEM_LOG("ATOM_TEST", Sev::inform)
506 << "Derivative hand calculation minus code " << nrm2_t_d_a;
507 if (nrm2_t_d_a > eps)
508 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
509 "This norm should be zero");
510 }
511
512 // Check second derivative
513 {
514 FTensor::Tensor2<double, 3, 3> t_S{
515
516 1., 1. / 2., 1. / 3.,
517
518 2. / 2., 1., 2. / 3.,
519
520 3. / 2., 1., 3. / 3.};
521
522 FTensor::Tensor2_symmetric<double, 3> t_S_sym;
523 t_S_sym(i, j) = t_S(i, j) || t_S(j, i);
524
525 auto t_dd =
526 EigenMatrix::getDiffDiffMat(t_L, t_N, f, d_f, dd_f, t_S_sym, 3);
527 auto t_dd_a = get_diff2_matrix2(t_S_sym, t_dd, FTensor::Number<3>());
528
529 MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_dd_a";
530 print_ddg(t_dd_a, "hand ");
531 MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_dd";
532 print_ddg(t_dd, "code ");
533
534 double nrm2_t_dd_a = get_norm_t4(t_dd_a, FTensor::Number<3>());
535 MOFEM_LOG("ATOM_TEST", Sev::inform)
536 << "Derivative hand calculation minus code " << nrm2_t_dd_a;
537 if (nrm2_t_dd_a > eps)
538 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
539 "This norm should be zero");
540 }
541 }
542 }
543
544 // Testing two equal eigenvalues
545 {
546
547 std::array<double, 9> a{5., 4., 0,
548
549 4., 5, 0.,
550
551 0.0, 0., 9};
552
553 auto tuple = run_lapack(a, swap01);
554 auto &t_a = std::get<0>(tuple);
555 auto &t_eig_vecs = std::get<1>(tuple);
556 auto &t_eig_vals = std::get<2>(tuple);
557
558 auto f = [](double v) { return v; };
559 auto d_f = [](double v) { return 1; };
560
561 constexpr double eps = 1e-10;
562 {
563 auto t_b = EigenMatrix::getMat(t_eig_vals, t_eig_vecs, f);
564 t_b(i, j) -= (t_a(i, j) || t_a(j, i)) / 2;
565 auto norm2_t_b = t_b(i, j) * t_b(i, j);
566 MOFEM_LOG("ATOM_TEST", Sev::inform)
567 << "Result should be matrix itself " << norm2_t_b;
568 if (norm2_t_b > eps)
569 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
570 "This norm should be zero");
571 }
572
573 {
574 auto t_d = EigenMatrix::getDiffMat(t_eig_vals, t_eig_vecs, f, d_f, 2);
575 auto t_d_a = get_diff_matrix(t_d, FTensor::Number<3>());
576 MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_d_a";
577 print_ddg(t_d_a, "hand ");
578 MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_d";
579 print_ddg(t_d, "code ");
580 double nrm2_t_d_a = get_norm_t4(t_d_a, FTensor::Number<3>());
581 MOFEM_LOG("ATOM_TEST", Sev::inform)
582 << "Derivative hand calculation minus code " << nrm2_t_d_a;
583 if (nrm2_t_d_a > eps)
584 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
585 "This norm should be zero");
586 }
587
588 {
589 auto f = [](double v) { return v * v; };
590 auto d_f = [](double v) { return 2 * v; };
591 auto t_d = EigenMatrix::getDiffMat(t_eig_vals, t_eig_vecs, f, d_f, 2);
592 auto t_d_a = get_diff_matrix2(t_a, t_d, FTensor::Number<3>());
593 MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_d_a";
594 print_ddg(t_d_a, "hand ");
595 MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_d";
596 print_ddg(t_d, "code ");
597 double nrm2_t_d_a = get_norm_t4(t_d_a, FTensor::Number<3>());
598 MOFEM_LOG("ATOM_TEST", Sev::inform)
599 << "Derivative hand calculation minus code " << nrm2_t_d_a;
600 if (nrm2_t_d_a > eps)
601 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
602 "This norm should be zero");
603 }
604 }
605
606 // Testing three equal eigenvalues
607 {
608
609 std::array<double, 9> a{4., 0., 0,
610
611 0., 4., 0.,
612
613 0.0, 0., 4.};
614
615 auto f = [](double v) { return v; };
616 auto d_f = [](double v) { return 1; };
617
618 auto tuple = run_lapack(a);
619 auto &t_a = std::get<0>(tuple);
620 auto &t_eig_vecs = std::get<1>(tuple);
621 auto &t_eig_vals = std::get<2>(tuple);
622
623 constexpr double eps = 1e-10;
624 {
625 auto t_b = EigenMatrix::getMat(t_eig_vals, t_eig_vecs, f);
626 t_b(i, j) -= (t_a(i, j) || t_a(j, i)) / 2;
627 auto norm2_t_b = t_b(i, j) * t_b(i, j);
628 MOFEM_LOG("ATOM_TEST", Sev::inform)
629 << "Result should be matrix itself " << norm2_t_b;
630 if (norm2_t_b > eps)
631 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
632 "This norm should be zero");
633 }
634
635 {
636 auto t_d = EigenMatrix::getDiffMat(t_eig_vals, t_eig_vecs, f, d_f, 1);
637 auto t_d_a = get_diff_matrix(t_d, FTensor::Number<3>());
638 MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_d_a";
639 print_ddg(t_d_a, "hand ");
640 MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_d";
641 print_ddg(t_d, "code ");
642 double nrm2_t_d_a = get_norm_t4(t_d_a, FTensor::Number<3>());
643 MOFEM_LOG("ATOM_TEST", Sev::inform)
644 << "Derivative hand calculation minus code " << nrm2_t_d_a;
645 if (nrm2_t_d_a > eps)
646 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
647 "This norm should be zero");
648 }
649
650 {
651 auto f = [](double v) { return v * v; };
652 auto d_f = [](double v) { return 2 * v; };
653 auto t_d = EigenMatrix::getDiffMat(t_eig_vals, t_eig_vecs, f, d_f, 1);
654 auto t_d_a = get_diff_matrix2(t_a, t_d, FTensor::Number<3>());
655 MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_d_a";
656 print_ddg(t_d_a, "hand ");
657 MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_d";
658 print_ddg(t_d, "code ");
659 double nrm2_t_d_a = get_norm_t4(t_d_a, FTensor::Number<3>());
660 MOFEM_LOG("ATOM_TEST", Sev::inform)
661 << "Derivative hand calculation minus code " << nrm2_t_d_a;
662 if (nrm2_t_d_a > eps)
663 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
664 "This norm should be zero");
665 }
666 }
667
668 // Check the second-derivative norm for a linear function near repeated
669 // eigenvalues
670 {
671
672 std::array<double, 9> a{0.1, 0., 0.,
673
674 0., 0.1, 0.,
675
676 0., 0., 0.1};
677
678 auto tuple = run_lapack(a);
679 // auto &t_a = std::get<0>(tuple);
680 auto &t_eig_vecs = std::get<1>(tuple);
681 auto &t_eig_vals = std::get<2>(tuple);
682
683 t_eig_vals(0) -= 1e-5;
684 t_eig_vals(2) += 1e-5;
685
686 constexpr double eps = 1e-10;
687
688 auto f = [](double v) { return v; };
689 auto d_f = [](double v) { return 1; };
690 auto dd_f = [](double v) { return 0; };
691
692 FTensor::Tensor2<double, 3, 3> t_S{
693
694 1., 1. / 2., 1. / 3.,
695
696 2. / 2., 1., 2. / 3.,
697
698 3. / 2., 1., 3. / 3.};
699
700 FTensor::Tensor2_symmetric<double, 3> t_S_sym;
701 t_S_sym(i, j) = t_S(i, j) || t_S(j, i);
702
703 auto t_dd = EigenMatrix::getDiffDiffMat(t_eig_vals, t_eig_vecs, f, d_f,
704 dd_f, t_S_sym, 1);
705
706 MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_dd";
707 print_ddg(t_dd, "test ");
708
709 double nrm2_t_dd = get_norm_t4(t_dd, FTensor::Number<3>());
710 MOFEM_LOG("ATOM_TEST", Sev::inform)
711 << "Second-derivative norm " << nrm2_t_dd;
712 if (nrm2_t_dd > eps)
713 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
714 "This norm should be zero");
715 }
716
717 // Check second derivative for a quadratic function with repeated
718 // eigenvalues
719 {
720
721 std::array<double, 9> a{2, 0., 0.,
722
723 0., 2, 0.,
724
725 0., 0., 2};
726
727 auto tuple = run_lapack(a);
728 // auto &t_a = std::get<0>(tuple);
729 auto &t_eig_vecs = std::get<1>(tuple);
730 auto &t_eig_vals = std::get<2>(tuple);
731
732 constexpr double eps = 1e-10;
733
734 auto f = [](double v) { return v * v; };
735 auto d_f = [](double v) { return 2 * v; };
736 auto dd_f = [](double v) { return 2; };
737 FTensor::Tensor2<double, 3, 3> t_S{
738
739 1., 1. / 2., 1. / 3.,
740
741 2. / 1., 1., 2. / 3.,
742
743 3. / 1., 3. / 1., 1.};
744
745 FTensor::Tensor2_symmetric<double, 3> t_S_sym;
746 t_S_sym(i, j) = t_S(i, j) || t_S(j, i);
747
748 auto t_dd = EigenMatrix::getDiffDiffMat(t_eig_vals, t_eig_vecs, f, d_f,
749 dd_f, t_S_sym, 1);
750 // print_ddg(t_dd, "test ");
751
752 auto t_dd_a = get_diff2_matrix2(t_S_sym, t_dd, FTensor::Number<3>());
753
754 double nrm2_t_dd_a = get_norm_t4(t_dd_a, FTensor::Number<3>());
755 MOFEM_LOG("ATOM_TEST", Sev::inform)
756 << "Derivative hand calculation minus code " << nrm2_t_dd_a;
757 if (nrm2_t_dd_a > eps)
758 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
759 "This norm should be zero");
760 }
761
762 // Check second derivative for two repeated eigenvalues
763 {
764
765 std::array<double, 9> a{5., 4., 0.,
766
767 4., 5., 0.,
768
769 0., 0., 9};
770
771 auto tuple = run_lapack(a, swap01);
772 // auto &t_a = std::get<0>(tuple);
773 auto &t_eig_vecs = std::get<1>(tuple);
774 auto &t_eig_vals = std::get<2>(tuple);
775
776 constexpr double eps = 1e-10;
777
778 auto f = [](double v) { return v * v; };
779 auto d_f = [](double v) { return 2 * v; };
780 auto dd_f = [](double v) { return 2; };
781
782 FTensor::Tensor2<double, 3, 3> t_S{
783
784 1., 1. / 2., 1. / 3.,
785
786 2. / 1., 1., 2. / 3.,
787
788 3. / 1., 3. / 1., 1.};
789
790 FTensor::Tensor2_symmetric<double, 3> t_S_sym;
791 t_S_sym(i, j) = t_S(i, j) || t_S(j, i);
792
793 auto t_dd = EigenMatrix::getDiffDiffMat(t_eig_vals, t_eig_vecs, f, d_f,
794 dd_f, t_S_sym, 2);
795 print_ddg(t_dd, "test ");
796
797 auto t_dd_a = get_diff2_matrix2(t_S_sym, t_dd, FTensor::Number<3>());
798
799 double nrm2_t_dd_a = get_norm_t4(t_dd_a, FTensor::Number<3>());
800 MOFEM_LOG("ATOM_TEST", Sev::inform)
801 << "Derivative hand calculation minus code " << nrm2_t_dd_a;
802 if (nrm2_t_dd_a > eps)
803 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
804 "This norm should be zero");
805 }
806
807 // Check second derivative for exponent near repeated eigenvalues
808 {
809
810 std::array<double, 9> a{2, 0., 0.,
811
812 0., 2, 0.,
813
814 0., 0., 2};
815
816 auto tuple = run_lapack(a);
817 // auto &t_a = std::get<0>(tuple);
818 auto &t_eig_vecs = std::get<1>(tuple);
819 auto &t_eig_vals = std::get<2>(tuple);
820
821 t_eig_vals(0) -= 1e-5;
822 t_eig_vals(2) += 1e-5;
823
824 constexpr double eps = 1e-7;
825
826 auto f = [](double v) { return exp(v); };
827 auto d_f = [](double v) { return exp(v); };
828 auto dd_f = [](double v) { return exp(v); };
829 FTensor::Tensor2<double, 3, 3> t_S{
830
831 1., 1. / 2., 1. / 3.,
832
833 2. / 1., 1., 2. / 3.,
834
835 3. / 1., 3. / 1., 1.};
836
837 FTensor::Tensor2_symmetric<double, 3> t_S_sym;
838 t_S_sym(i, j) = t_S(i, j) || t_S(j, i);
839
840 auto t_dd_1 = EigenMatrix::getDiffDiffMat(t_eig_vals, t_eig_vecs, f, d_f,
841 dd_f, t_S_sym, 3);
842 auto t_dd_2 = EigenMatrix::getDiffDiffMat(t_eig_vals, t_eig_vecs, f, d_f,
843 dd_f, t_S_sym, 1);
844
845 double nrm2_t_dd_t1 = get_norm_t4(t_dd_1, FTensor::Number<3>());
846 MOFEM_LOG("ATOM_TEST", Sev::verbose)
847 << "Derivative norm t_dd_1 " << nrm2_t_dd_t1;
848
849 double nrm2_t_dd_t2 = get_norm_t4(t_dd_2, FTensor::Number<3>());
850 MOFEM_LOG("ATOM_TEST", Sev::verbose)
851 << "Derivative norm t_dd_2 " << nrm2_t_dd_t2;
852
853 print_ddg(t_dd_1, "t_dd_1 ");
854 print_ddg(t_dd_2, "t_dd_2 ");
855
856 FTensor::Ddg<double, 3, 3> t_dd_3;
857 t_dd_3(i, j, k, l) = t_dd_1(i, j, k, l) - t_dd_2(i, j, k, l);
858
859 for (int ii = 0; ii != 3; ++ii)
860 for (int jj = 0; jj != 3; ++jj)
861 for (int kk = 0; kk != 3; ++kk)
862 for (int ll = 0; ll != 3; ++ll) {
863 constexpr double eps = 1e-4;
864 if (std::abs(t_dd_3(ii, jj, kk, ll)) > eps)
865 MOFEM_LOG("ATOM_TEST", Sev::error)
866 << "Error " << ii << " " << jj << " " << kk << " " << ll
867 << " " << t_dd_1(ii, jj, kk, ll) << " "
868 << t_dd_2(ii, jj, kk, ll);
869 }
870
871 double nrm2_t_dd_3 = get_norm_t4(t_dd_3, FTensor::Number<3>());
872 MOFEM_LOG("ATOM_TEST", Sev::inform)
873 << "Derivative branch difference " << nrm2_t_dd_3;
874 if (nrm2_t_dd_3 > eps)
875 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
876 "This norm should be zero");
877 }
878
879 // Check second derivative for quadratic function against perturbed
880 // repeated eigenvalues
881 {
882
883 std::array<double, 9> a{5., 4., 0.,
884
885 4., 5., 0.,
886
887 0., 0., 9};
888
889 auto tuple = run_lapack(a, swap01);
890 // auto &t_a = std::get<0>(tuple);
891 auto &t_eig_vecs = std::get<1>(tuple);
892 auto &t_eig_vals = std::get<2>(tuple);
893
894 t_eig_vals(0) -= 1e-4;
895 t_eig_vals(2) += 1e-4;
896
897 constexpr double eps = 1e-4;
898
899 auto f = [](double v) { return v * v; };
900 auto d_f = [](double v) { return 2 * v; };
901 auto dd_f = [](double v) { return 2; };
902
903 FTensor::Tensor2<double, 3, 3> t_S{
904
905 1., 1. / 2., 1. / 3.,
906
907 2. / 1., 1., 2. / 3.,
908
909 3. / 1., 3. / 1., 1.};
910
911 FTensor::Tensor2_symmetric<double, 3> t_S_sym;
912 t_S_sym(i, j) = t_S(i, j) || t_S(j, i);
913
914 auto t_dd_1 = EigenMatrix::getDiffDiffMat(t_eig_vals, t_eig_vecs, f, d_f,
915 dd_f, t_S_sym, 3);
916 auto t_dd_2 = EigenMatrix::getDiffDiffMat(t_eig_vals, t_eig_vecs, f, d_f,
917 dd_f, t_S_sym, 2);
918
919 double nrm2_t_dd_t1 = get_norm_t4(t_dd_1, FTensor::Number<3>());
920 MOFEM_LOG("ATOM_TEST", Sev::verbose)
921 << "Derivative norm t_dd_1 " << nrm2_t_dd_t1;
922
923 double nrm2_t_dd_t2 = get_norm_t4(t_dd_2, FTensor::Number<3>());
924 MOFEM_LOG("ATOM_TEST", Sev::verbose)
925 << "Derivative norm t_dd_2 " << nrm2_t_dd_t2;
926
927 print_ddg(t_dd_1, "t_dd_1 ");
928 print_ddg(t_dd_2, "t_dd_2 ");
929
930 FTensor::Ddg<double, 3, 3> t_dd_3;
931 t_dd_3(i, j, k, l) = t_dd_1(i, j, k, l) - t_dd_2(i, j, k, l);
932
933 for (int ii = 0; ii != 3; ++ii)
934 for (int jj = 0; jj != 3; ++jj)
935 for (int kk = 0; kk != 3; ++kk)
936 for (int ll = 0; ll != 3; ++ll) {
937 constexpr double eps = 1e-3;
938 if (std::abs(t_dd_3(ii, jj, kk, ll)) > eps)
939 MOFEM_LOG("ATOM_TEST", Sev::error)
940 << "Error " << ii << " " << jj << " " << kk << " " << ll
941 << " " << t_dd_1(ii, jj, kk, ll) << " "
942 << t_dd_2(ii, jj, kk, ll);
943 }
944
945 double nrm2_t_dd_3 = get_norm_t4(t_dd_3, FTensor::Number<3>());
946 MOFEM_LOG("ATOM_TEST", Sev::inform)
947 << "Derivative branch difference " << nrm2_t_dd_3;
948 if (nrm2_t_dd_3 > eps)
949 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
950 "This norm should be zero");
951 }
952
953 // Check second derivative for a cubic polynomial near repeated
954 // eigenvalues
955 {
956
957 std::array<double, 9> a{5., 4., 0.,
958
959 4., 5., 0.,
960
961 0., 0., 9};
962
963 auto tuple = run_lapack(a, swap01);
964 // auto &t_a = std::get<0>(tuple);
965 auto &t_eig_vecs = std::get<1>(tuple);
966 auto &t_eig_vals = std::get<2>(tuple);
967
968 constexpr double eps = 1e-4;
969 constexpr int p = 3;
970
971 auto f = [](double v) { return pow(v, p); };
972 auto d_f = [](double v) { return p * pow(v, p - 1); };
973 auto dd_f = [](double v) {
974 return p * (p - 1) * pow(v, std::max(0, p - 2));
975 };
976
977 FTensor::Tensor2<double, 3, 3> t_S{
978
979 1., 1. / 2., 1. / 3.,
980
981 2. / 1., 1., 2. / 3.,
982
983 3. / 1., 3. / 1., 1.};
984
985 // FTensor::Tensor2_symmetric<double, 3> t_S_sym;
986 // t_S_sym(i, j) = t_S(i, j) || t_S(j, i);
987
988 t_eig_vals(0) += 2e-5;
989 t_eig_vals(2) -= 2e-5;
990 auto t_dd_1 = EigenMatrix::getDiffDiffMat(t_eig_vals, t_eig_vecs, f, d_f,
991 dd_f, t_S, 3);
992 auto t_dd_2 = EigenMatrix::getDiffDiffMat(t_eig_vals, t_eig_vecs, f, d_f,
993 dd_f, t_S, 2);
994
995 double nrm2_t_dd_t1 = get_norm_t4(t_dd_1, FTensor::Number<3>());
996 MOFEM_LOG("ATOM_TEST", Sev::verbose)
997 << "Derivative norm t_dd_1 " << nrm2_t_dd_t1;
998
999 double nrm2_t_dd_t2 = get_norm_t4(t_dd_2, FTensor::Number<3>());
1000 MOFEM_LOG("ATOM_TEST", Sev::verbose)
1001 << "Derivative norm t_dd_2 " << nrm2_t_dd_t2;
1002
1003 print_ddg(t_dd_1, "t_dd_1 ");
1004 print_ddg(t_dd_2, "t_dd_2 ");
1005
1006 FTensor::Ddg<double, 3, 3> t_dd_3;
1007 t_dd_3(i, j, k, l) = t_dd_1(i, j, k, l) - t_dd_2(i, j, k, l);
1008
1009 for (int ii = 0; ii != 3; ++ii)
1010 for (int jj = 0; jj != 3; ++jj)
1011 for (int kk = 0; kk != 3; ++kk)
1012 for (int ll = 0; ll != 3; ++ll) {
1013 constexpr double eps = 1e-3;
1014 if (std::abs(t_dd_3(ii, jj, kk, ll)) > eps)
1015 MOFEM_LOG("ATOM_TEST", Sev::error)
1016 << "Error " << ii << " " << jj << " " << kk << " " << ll
1017 << " " << t_dd_1(ii, jj, kk, ll) << " "
1018 << t_dd_2(ii, jj, kk, ll) << " " << t_dd_3(ii, jj, kk, ll);
1019 }
1020
1021 double nrm2_t_dd_3 = get_norm_t4(t_dd_3, FTensor::Number<3>());
1022 MOFEM_LOG("ATOM_TEST", Sev::inform)
1023 << "Derivative branch difference " << nrm2_t_dd_3;
1024 if (nrm2_t_dd_3 > eps)
1025 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1026 "This norm should be zero");
1027 }
1028
1029 // Speed
1030 {
1031
1032 std::array<double, 9> a{1., 0.1, -0.5,
1033
1034 0.1, 2., 0.,
1035
1036 -0.5, 0., 3.};
1037
1038 auto tuple = run_lapack(a);
1039 // auto &t_a = std::get<0>(tuple);
1040 auto &t_eig_vecs = std::get<1>(tuple);
1041 auto &t_eig_vals = std::get<2>(tuple);
1042
1043 auto f = [](double v) { return exp(v); };
1044 auto d_f = [](double v) { return exp(v); };
1045 auto dd_f = [](double v) { return exp(v); };
1046
1047 FTensor::Tensor2<double, 3, 3> t_S{
1048
1049 1., 1. / 2., 1. / 3.,
1050
1051 2. / 1., 1., 2. / 3.,
1052
1053 3. / 1., 3. / 1., 1.};
1054
1055 FTensor::Tensor2_symmetric<double, 3> t_S_sym;
1056 t_S_sym(i, j) = t_S(i, j) || t_S(j, i);
1057
1058 MOFEM_LOG("ATOM_TEST", Sev::inform) << "Start";
1059 for (int ii = 0; ii != 1000; ++ii) {
1060 auto t_d = EigenMatrix::getDiffMat(t_eig_vals, t_eig_vecs, f, d_f, 3);
1061 auto t_dd = EigenMatrix::getDiffDiffMat(t_eig_vals, t_eig_vecs, f, d_f,
1062 dd_f, t_S_sym, 3);
1063 std::ignore = t_d;
1064 std::ignore = t_dd;
1065 }
1066 MOFEM_LOG("ATOM_TEST", Sev::inform) << "End";
1067 }
1068
1069 // 2d case
1070
1071 auto run_lapack_2d = [](auto &a) {
1072 int info;
1073 double wkopt;
1074 double w[2];
1075
1076 FTensor::Tensor2<double, 2, 2> t_a{
1077
1078 a[0], a[1],
1079
1080 a[2], a[3]};
1081
1082 /* Query and allocate the optimal workspace */
1083 int lwork = -1;
1084 info = lapack_dsyev('V', 'U', 2, a.data(), 2, w, &wkopt, lwork);
1085 if (info > 0)
1086 THROW_MESSAGE("The algorithm failed to compute eigenvalues.");
1087 lwork = (int)wkopt;
1088 std::vector<double> work(lwork);
1089 /* Solve eigenproblem */
1090 info = lapack_dsyev('V', 'U', 2, a.data(), 2, w, &*work.begin(), lwork);
1091 if (info > 0)
1092 THROW_MESSAGE("The algorithm failed to compute eigenvalues.");
1093
1094 FTensor::Tensor2<double, 2, 2> t_eig_vecs{
1095
1096 a[0 * 2 + 0], a[0 * 2 + 1],
1097
1098 a[1 * 2 + 0], a[1 * 2 + 1]};
1099
1100 FTensor::Tensor1<double, 2> t_eig_vals{w[0], w[1]};
1101
1102 return std::make_tuple(t_a, t_eig_vecs, t_eig_vals);
1103 };
1104
1105 // Testing quadratic function for 2d
1106 {
1107
1108 std::array<double, 9> a{1., 0.1,
1109
1110 0.1, 2.};
1111
1112 auto tuple = run_lapack_2d(a);
1113 auto &t_A = std::get<0>(tuple);
1114 auto &t_eig_vecs = std::get<1>(tuple);
1115 auto &t_eig_vals = std::get<2>(tuple);
1116
1117 auto f = [](double v) { return v * v; };
1118 auto d_f = [](double v) { return 2 * v; };
1119 auto dd_f = [](double v) { return 2; };
1120
1121 constexpr double eps = 1e-6;
1122
1123 FTensor::Index<'i', 2> i;
1124 FTensor::Index<'j', 2> j;
1125 FTensor::Index<'k', 2> k;
1126 FTensor::Index<'l', 2> l;
1127
1128 // check if multiplication gives right value
1129 {
1130 auto t_b = EigenMatrix::getMat(t_eig_vals, t_eig_vecs, f);
1131 FTensor::Tensor2<double, 2, 2> t_a;
1132 t_a(i, j) = t_b(i, j) - t_A(i, k) * t_A(k, j);
1133 print_mat(t_a);
1134 auto norm2_t_a = t_a(i, j) * t_a(i, j);
1135 MOFEM_LOG("ATOM_TEST", Sev::inform)
1136 << "Result should be matrix times matrix " << norm2_t_a;
1137 if (norm2_t_a > eps)
1138 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1139 "This norm should be zero");
1140 }
1141
1142 // Check first derivative
1143 {
1144 auto t_d = EigenMatrix::getDiffMat(t_eig_vals, t_eig_vecs, f, d_f, 2);
1145 auto t_d_a = get_diff_matrix2(t_A, t_d, FTensor::Number<2>());
1146 double nrm2_t_d_a = get_norm_t4(t_d_a, FTensor::Number<2>());
1147 MOFEM_LOG("ATOM_TEST", Sev::inform)
1148 << "Derivative hand calculation minus code " << nrm2_t_d_a;
1149 if (nrm2_t_d_a > eps)
1150 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1151 "This norm should be zero");
1152 }
1153
1154 // Check second derivative
1155 {
1156 FTensor::Tensor2<double, 2, 2> t_S{
1157
1158 1., 1. / 2,
1159
1160 2. / 2., 1.};
1161
1162 // FTensor::Index<'i', 2> i;
1163 // FTensor::Index<'j', 2> j;
1164 // FTensor::Tensor2_symmetric<double, 2> t_S_sym;
1165 // t_S_sym(i, j) = t_S(i, j) || t_S(j, i);
1166
1167 auto t_dd = EigenMatrix::getDiffDiffMat(t_eig_vals, t_eig_vecs, f, d_f,
1168 dd_f, t_S, 2);
1169 auto t_dd_a = get_diff2_matrix2(t_S, t_dd, FTensor::Number<2>());
1170
1171 MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_dd_a";
1172 MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_dd";
1173
1174 double nrm2_t_dd_a = get_norm_t4(t_dd_a, FTensor::Number<2>());
1175 MOFEM_LOG("ATOM_TEST", Sev::inform)
1176 << "Derivative hand calculation minus code " << nrm2_t_dd_a;
1177 if (nrm2_t_dd_a > eps)
1178 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1179 "This norm should be zero");
1180 }
1181 }
1182
1183 // Testing quadratic function for repeating eigenvalues
1184 {
1185
1186 std::array<double, 9> a{2., 0,
1187
1188 0, 2.};
1189
1190 auto tuple = run_lapack_2d(a);
1191 auto &t_A = std::get<0>(tuple);
1192 auto &t_eig_vecs = std::get<1>(tuple);
1193 auto &t_eig_vals = std::get<2>(tuple);
1194
1195 auto f = [](double v) { return v * v; };
1196 auto d_f = [](double v) { return 2 * v; };
1197 auto dd_f = [](double v) { return 2; };
1198
1199 constexpr double eps = 1e-6;
1200
1201 FTensor::Index<'i', 2> i;
1202 FTensor::Index<'j', 2> j;
1203 FTensor::Index<'k', 2> k;
1204 FTensor::Index<'l', 2> l;
1205
1206 // check if multiplication gives right value
1207 {
1208 auto t_b = EigenMatrix::getMat(t_eig_vals, t_eig_vecs, f);
1209 FTensor::Tensor2<double, 2, 2> t_a;
1210 t_a(i, j) = t_b(i, j) - t_A(i, k) * t_A(k, j);
1211 print_mat(t_a);
1212 auto norm2_t_a = t_a(i, j) * t_a(i, j);
1213 MOFEM_LOG("ATOM_TEST", Sev::inform)
1214 << "Result should be matrix times matrix " << norm2_t_a;
1215 if (norm2_t_a > eps)
1216 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1217 "This norm should be zero");
1218 }
1219
1220 // Check first derivative
1221 {
1222 auto t_d = EigenMatrix::getDiffMat(t_eig_vals, t_eig_vecs, f, d_f, 1);
1223 auto t_d_a = get_diff_matrix2(t_A, t_d, FTensor::Number<2>());
1224 double nrm2_t_d_a = get_norm_t4(t_d_a, FTensor::Number<2>());
1225 MOFEM_LOG("ATOM_TEST", Sev::inform)
1226 << "Derivative hand calculation minus code " << nrm2_t_d_a;
1227 if (nrm2_t_d_a > eps)
1228 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1229 "This norm should be zero");
1230 }
1231
1232 // Check second derivative
1233 {
1234 FTensor::Tensor2<double, 2, 2> t_S{
1235
1236 1., 1. / 2,
1237
1238 2. / 2., 1.};
1239
1240 FTensor::Index<'i', 2> i;
1241 FTensor::Index<'j', 2> j;
1242 FTensor::Tensor2_symmetric<double, 2> t_S_sym;
1243 t_S_sym(i, j) = t_S(i, j) || t_S(j, i);
1244
1245 auto t_dd = EigenMatrix::getDiffDiffMat(t_eig_vals, t_eig_vecs, f, d_f,
1246 dd_f, t_S_sym, 1);
1247 auto t_dd_a = get_diff2_matrix2(t_S_sym, t_dd, FTensor::Number<2>());
1248
1249 MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_dd_a";
1250 MOFEM_LOG("ATOM_TEST", Sev::verbose) << "t_dd";
1251
1252 double nrm2_t_dd_a = get_norm_t4(t_dd_a, FTensor::Number<2>());
1253 MOFEM_LOG("ATOM_TEST", Sev::inform)
1254 << "Derivative hand calculation minus code " << nrm2_t_dd_a;
1255 if (nrm2_t_dd_a > eps)
1256 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1257 "This norm should be zero");
1258 }
1259 }
1260 }
1261 CATCH_ERRORS;
1262
1263 CHKERR MoFEM::Core::Finalize();
1264}
FTensor.hpp
Tensors class implemented by Walter Landry.
main
int main()
Definition adol-c_atom.cpp:46
a
constexpr double a
Definition approx_sphere.cpp:30
eps
static const double eps
Definition check_base_functions_derivatives_on_tet.cpp:11
FTensor::Kronecker_Delta_symmetric
Kronecker Delta class symmetric.
Definition Kronecker_Delta.hpp:49
FTensor::Kronecker_Delta
Kronecker Delta class.
Definition Kronecker_Delta.hpp:15
FTensor::Tensor2_symmetric
Definition Tensor2_symmetric_value.hpp:14
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition definitions.h:385
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition definitions.h:40
CHKERR
#define CHKERR
Inline error check.
Definition definitions.h:551
THROW_MESSAGE
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
Definition definitions.h:577
print_mat
void print_mat(double *M, int m, int n)
print matric M
t_kd
constexpr auto t_kd
Definition free_surface.cpp:161
MoFEM::LogManager::setLog
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition LogManager.cpp:392
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition LogManager.hpp:316
MOFEM_LOG_ATTRIBUTES
#define MOFEM_LOG_ATTRIBUTES(channel, bit)
Add attributes to channel.
Definition LogManager.hpp:304
MoFEM::LogManager::BitLineID
@ BitLineID
Definition LogManager.hpp:48
v
const double v
phase velocity of light in medium (cm/ns)
Definition initial_diffusion.cpp:40
lapack_dsyev
static __CLPK_integer lapack_dsyev(char jobz, char uplo, __CLPK_integer n, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_doublereal *w, __CLPK_doublereal *work, __CLPK_integer lwork)
Definition lapack_wrap.h:273
get_diff_matrix
auto get_diff_matrix(T1 &t_d, const FTensor::Number< DIM > &)
Definition matrix_function.cpp:142
diff_ddg
void diff_ddg(T1 &t_1, T2 &t_2, const FTensor::Number< DIM > &)
Definition matrix_function.cpp:26
help
static char help[]
Definition matrix_function.cpp:182
l
FTensor::Index< 'l', 3 > l
Definition matrix_function.cpp:23
j
FTensor::Index< 'j', 3 > j
Definition matrix_function.cpp:21
get_diff2_matrix2
auto get_diff2_matrix2(T1 &t_s, T2 &t_dd, const FTensor::Number< DIM > &)
Definition matrix_function.cpp:84
i
FTensor::Index< 'i', 3 > i
Definition matrix_function.cpp:20
get_norm_t4
auto get_norm_t4(T &t, const FTensor::Number< DIM > &)
Definition matrix_function.cpp:172
k
FTensor::Index< 'k', 3 > k
Definition matrix_function.cpp:22
get_diff_matrix2
auto get_diff_matrix2(T1 &t_a, T2 &t_d, const FTensor::Number< DIM > &)
Definition matrix_function.cpp:47
EigenMatrix::getMat
auto getMat(A &&t_val, B &&t_vec, Fun< double > f)
Get the Mat object.
Definition MatrixFunction.hpp:133
EigenMatrix::getDiffMat
auto getDiffMat(A &&t_val, B &&t_vec, Fun< double > f, Fun< double > d_f, const int nb)
Get the Diff Mat object.
Definition MatrixFunction.hpp:201
EigenMatrix::getDiffDiffMat
auto getDiffDiffMat(A &&t_val, B &&t_vec, Fun< double > f, Fun< double > d_f, Fun< double > dd_f, C &&t_S, const int nb)
Get the Diff Diff Mat object.
Definition MatrixFunction.hpp:414
MoFEM
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
t
constexpr double t
plate stiffness
Definition plate.cpp:58
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition Core.cpp:68
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition Core.cpp:123
MoFEM::LogManager::createSink
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition LogManager.cpp:293
MoFEM::LogManager::getStrmSelf
static boost::shared_ptr< std::ostream > getStrmSelf()
Get the strm self object.
Definition LogManager.cpp:335
Generated by Doxygen 1.9.8 and hosted at University of Glasgow