--- a/experiments/tests/test_minimisation.py Mon Jun 29 00:20:42 2026 +0300
+++ b/experiments/tests/test_minimisation.py Sun Oct 04 13:02:26 2026 -0500
@@ -1,9 +1,23 @@
+import numpy as np
+from dolfinx import fem, mesh
+from mpi4py import MPI
+
from dolfinx_access import eval_Function_f64, min_Function_f64_p2
# A simple test that the finite-dimensional subproblem is solved correctly
-def test_minimisation(pde):
- test = fem.Function(pde.V)
+def test_minimisation(domain_dims=[[-5, 5], [-5, 5]], degree=2, nx=256, ny=256, n=512):
+ domain = mesh.create_rectangle(
+ MPI.COMM_SELF,
+ [
+ np.array([domain_dims[0][0], domain_dims[1][0]]),
+ np.array([domain_dims[0][1], domain_dims[1][1]]),
+ ],
+ [nx, ny],
+ mesh.CellType.triangle,
+ )
+ V = fem.functionspace(domain, ("Lagrange", degree))
+ test = fem.Function(V)
# Beale function
def beale(x):
@@ -13,22 +27,69 @@
+ (2.625 - x[0] + x[0] * (x[1] ** 3)) ** 2
)
- test.interpolate(beale)
- m = min_Function_f64_p2(test._cpp_object)
- print("Beale minimizer:", m.x, m.v)
-
- def quadratic_testf(x):
+ def quadratic(x):
return x[0] ** 2 - 5 * x[1] ** 2 + x[0] * x[1] - x[0] + x[1]
- testf = quadratic_testf
- test.interpolate(testf)
- n = 256
- dist = 0.0
- for x in np.linspace(pde.domain_size[0], pde.domain_size[1], n):
- for y in np.linspace(pde.domain_size[2], pde.domain_size[3], n):
- p = [x, y, 0.0]
- v = eval_Function_f64(test._cpp_object, p)
- v0 = testf(p)
- d = np.linalg.norm(v - v0)
- dist = max(d, dist)
- print("Quadratic function evaluation max discrepancy:", dist)
+ def sphere(x):
+ return x[0] ** 2 + x[1] ** 2
+
+ def booth(x):
+ return (x[0] + 2 * x[1] - 7) ** 2 + (2 * x[0] + x[1] - 5) ** 2
+
+ def himmelblau(x):
+ return (x[0] ** 2 + x[1] - 11) ** 2 + (x[0] + x[1] ** 2 - 7) ** 2
+
+ for testf, name in [
+ (beale, "Beale"),
+ (quadratic, "Quadratic"),
+ (sphere, "Sphere"),
+ (booth, "Booth"),
+ (himmelblau, "Himmelblau"),
+ ]:
+ test.interpolate(testf)
+ m = min_Function_f64_p2(test._cpp_object)
+ print(
+ name,
+ "minimizer:",
+ m.v,
+ "(verification ",
+ eval_Function_f64(test._cpp_object, [m.x[0], m.x[1], 0]),
+ " original ",
+ testf(m.x),
+ ") @ ",
+ m.x,
+ )
+
+ dist = 0.0
+ m = float("inf")
+ m0 = float("inf")
+ mp = None
+ for x in np.linspace(domain_dims[0][0], domain_dims[0][1], n):
+ for y in np.linspace(domain_dims[1][0], domain_dims[1][1], n):
+ p = [x, y, 0.0]
+ v = eval_Function_f64(test._cpp_object, p)
+ v0 = testf(p)
+ if v0 < m0:
+ m = v
+ m0 = v0
+ mp = p
+ d = np.linalg.norm(v - v0)
+ dist = max(d, dist)
+ print(name, "evaluation max discrepancy:", dist)
+ print(
+ name,
+ "rough scan minimiser ",
+ m,
+ "(model",
+ m0,
+ ") @ [",
+ float(mp[0]),
+ ", ",
+ float(mp[1]),
+ "]",
+ )
+ print()
+
+
+test_minimisation()
+exit(0)