Mercurial > repos > pointsource_pde / file comparison
comparison: experiments/tests/test_minimisation.py
experiments/tests/test_minimisation.py
- changeset 12
- a0d61ac3719e
- parent 1
- a4137aedcb3a
equal
deleted
inserted
replaced
| 11:c4f923017349 | 12:a0d61ac3719e |
|---|---|
| 1 import numpy as np | |
| 2 from dolfinx import fem, mesh | |
| 3 from mpi4py import MPI | |
| 4 | |
| 1 from dolfinx_access import eval_Function_f64, min_Function_f64_p2 | 5 from dolfinx_access import eval_Function_f64, min_Function_f64_p2 |
| 2 | 6 |
| 3 | 7 |
| 4 # A simple test that the finite-dimensional subproblem is solved correctly | 8 # A simple test that the finite-dimensional subproblem is solved correctly |
| 5 def test_minimisation(pde): | 9 def test_minimisation(domain_dims=[[-5, 5], [-5, 5]], degree=2, nx=256, ny=256, n=512): |
| 6 test = fem.Function(pde.V) | 10 domain = mesh.create_rectangle( |
| 11 MPI.COMM_SELF, | |
| 12 [ | |
| 13 np.array([domain_dims[0][0], domain_dims[1][0]]), | |
| 14 np.array([domain_dims[0][1], domain_dims[1][1]]), | |
| 15 ], | |
| 16 [nx, ny], | |
| 17 mesh.CellType.triangle, | |
| 18 ) | |
| 19 V = fem.functionspace(domain, ("Lagrange", degree)) | |
| 20 test = fem.Function(V) | |
| 7 | 21 |
| 8 # Beale function | 22 # Beale function |
| 9 def beale(x): | 23 def beale(x): |
| 10 return ( | 24 return ( |
| 11 (1.5 - x[0] + x[0] * x[1]) ** 2 | 25 (1.5 - x[0] + x[0] * x[1]) ** 2 |
| 12 + (2.25 - x[0] + x[0] * (x[1] ** 2)) ** 2 | 26 + (2.25 - x[0] + x[0] * (x[1] ** 2)) ** 2 |
| 13 + (2.625 - x[0] + x[0] * (x[1] ** 3)) ** 2 | 27 + (2.625 - x[0] + x[0] * (x[1] ** 3)) ** 2 |
| 14 ) | 28 ) |
| 15 | 29 |
| 16 test.interpolate(beale) | 30 def quadratic(x): |
| 17 m = min_Function_f64_p2(test._cpp_object) | |
| 18 print("Beale minimizer:", m.x, m.v) | |
| 19 | |
| 20 def quadratic_testf(x): | |
| 21 return x[0] ** 2 - 5 * x[1] ** 2 + x[0] * x[1] - x[0] + x[1] | 31 return x[0] ** 2 - 5 * x[1] ** 2 + x[0] * x[1] - x[0] + x[1] |
| 22 | 32 |
| 23 testf = quadratic_testf | 33 def sphere(x): |
| 24 test.interpolate(testf) | 34 return x[0] ** 2 + x[1] ** 2 |
| 25 n = 256 | 35 |
| 26 dist = 0.0 | 36 def booth(x): |
| 27 for x in np.linspace(pde.domain_size[0], pde.domain_size[1], n): | 37 return (x[0] + 2 * x[1] - 7) ** 2 + (2 * x[0] + x[1] - 5) ** 2 |
| 28 for y in np.linspace(pde.domain_size[2], pde.domain_size[3], n): | 38 |
| 29 p = [x, y, 0.0] | 39 def himmelblau(x): |
| 30 v = eval_Function_f64(test._cpp_object, p) | 40 return (x[0] ** 2 + x[1] - 11) ** 2 + (x[0] + x[1] ** 2 - 7) ** 2 |
| 31 v0 = testf(p) | 41 |
| 32 d = np.linalg.norm(v - v0) | 42 for testf, name in [ |
| 33 dist = max(d, dist) | 43 (beale, "Beale"), |
| 34 print("Quadratic function evaluation max discrepancy:", dist) | 44 (quadratic, "Quadratic"), |
| 45 (sphere, "Sphere"), | |
| 46 (booth, "Booth"), | |
| 47 (himmelblau, "Himmelblau"), | |
| 48 ]: | |
| 49 test.interpolate(testf) | |
| 50 m = min_Function_f64_p2(test._cpp_object) | |
| 51 print( | |
| 52 name, | |
| 53 "minimizer:", | |
| 54 m.v, | |
| 55 "(verification ", | |
| 56 eval_Function_f64(test._cpp_object, [m.x[0], m.x[1], 0]), | |
| 57 " original ", | |
| 58 testf(m.x), | |
| 59 ") @ ", | |
| 60 m.x, | |
| 61 ) | |
| 62 | |
| 63 dist = 0.0 | |
| 64 m = float("inf") | |
| 65 m0 = float("inf") | |
| 66 mp = None | |
| 67 for x in np.linspace(domain_dims[0][0], domain_dims[0][1], n): | |
| 68 for y in np.linspace(domain_dims[1][0], domain_dims[1][1], n): | |
| 69 p = [x, y, 0.0] | |
| 70 v = eval_Function_f64(test._cpp_object, p) | |
| 71 v0 = testf(p) | |
| 72 if v0 < m0: | |
| 73 m = v | |
| 74 m0 = v0 | |
| 75 mp = p | |
| 76 d = np.linalg.norm(v - v0) | |
| 77 dist = max(d, dist) | |
| 78 print(name, "evaluation max discrepancy:", dist) | |
| 79 print( | |
| 80 name, | |
| 81 "rough scan minimiser ", | |
| 82 m, | |
| 83 "(model", | |
| 84 m0, | |
| 85 ") @ [", | |
| 86 float(mp[0]), | |
| 87 ", ", | |
| 88 float(mp[1]), | |
| 89 "]", | |
| 90 ) | |
| 91 print() | |
| 92 | |
| 93 | |
| 94 test_minimisation() | |
| 95 exit(0) |
