{ "cells": [ { "cell_type": "markdown", "id": "467ba05b", "metadata": {}, "source": [ "# Thickness optimization and shape derivatives\n", "\n", "Two things the analysis walkthrough only mentions: driving an optimizer with Hermit's\n", "outputs, and differentiating with respect to the **mesh coordinates**.\n", "\n", "Both run here for real, on a small cantilever plate. The optimization needs the\n", "optional solver stack (`pip install 'hermit[opt]'`)." ] }, { "cell_type": "markdown", "id": "3787ac35", "metadata": {}, "source": [ "## Setup\n", "\n", "The same 2 \u00d7 10 cantilever plate as the analysis walkthrough, coarser so the optimizer\n", "finishes quickly." ] }, { "cell_type": "code", "execution_count": 1, "id": "23b591dd", "metadata": { "execution": { "iopub.execute_input": "2026-09-09T21:40:21.082344Z", "iopub.status.busy": "2026-09-09T21:40:21.082256Z", "iopub.status.idle": "2026-09-09T21:40:21.750497Z", "shell.execute_reply": "2026-09-09T21:40:21.750018Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "33 vertices\n" ] } ], "source": [ "import inspect\n", "\n", "import numpy as np\n", "import basix.ufl\n", "import ufl\n", "from mpi4py import MPI\n", "from dolfinx.mesh import create_mesh\n", "import csdl_alpha as csdl\n", "\n", "import hermit as hm\n", "\n", "LENGTH, WIDTH = 10.0, 2.0\n", "E_VAL, NU_VAL, H_VAL, RHO_VAL = 4.32e8, 0.0, 0.2, 1.0\n", "PRESSURE_Z = 2.0\n", "nx, ny = 10, 2\n", "\n", "\n", "def plate_mesh(nx, ny):\n", " xs = np.linspace(0.0, LENGTH, nx + 1)\n", " ys = np.linspace(0.0, WIDTH, ny + 1)\n", " points = np.array([[x, y, 0.0] for y in ys for x in xs])\n", " cells = np.array([[j * (nx + 1) + i, j * (nx + 1) + i + 1,\n", " (j + 1) * (nx + 1) + i, (j + 1) * (nx + 1) + i + 1]\n", " for j in range(ny) for i in range(nx)], dtype=np.int64)\n", " el = ufl.Mesh(basix.ufl.element(\"Lagrange\", \"quadrilateral\", 1, shape=(3,)))\n", " if list(inspect.signature(create_mesh).parameters)[2] == \"e\": # DOLFINx 0.11\n", " return create_mesh(MPI.COMM_WORLD, cells, el, points)\n", " return create_mesh(MPI.COMM_WORLD, cells, points, el) # DOLFINx 0.9\n", "\n", "\n", "mesh = plate_mesh(nx, ny)\n", "print(f\"{mesh.geometry.x.shape[0]} vertices\")" ] }, { "cell_type": "markdown", "id": "afe5fb57", "metadata": {}, "source": [ "## Thickness optimization\n", "\n", "Minimise compliance at fixed mass. The whole model is one CSDL graph, so the optimizer\n", "interface is three calls on the variables \u2014 `set_as_design_variable`,\n", "`set_as_constraint`, `set_as_objective` \u2014 and then handing the recorder to modopt.\n", "\n", "**Choose the design space deliberately.** Giving every node its own free thickness is\n", "the obvious thing to do and it is ill-posed: bending stiffness goes as $\\int t^3$ while\n", "mass goes as $\\int t$, so at fixed mass an *oscillating* thickness genuinely raises the\n", "stiffness. A free nodal design space therefore has mesh-dependent checkerboard optima\n", "that really do beat any smooth profile in the discrete objective \u2014 the optimizer is\n", "right, the target is wrong. The usual remedies are a filter, a perimeter penalty, or a\n", "restricted design space.\n", "\n", "We take the last, and it makes the example sharper. Restrict the thickness to\n", "$t(x) = A\\,\\xi^{p}$ with $\\xi = (L-x)/L + \\varepsilon$, and let the optimizer find both\n", "$A$ and the exponent $p$. In the beam limit $C = \\int M^2/t^3\\,\\mathrm{d}x$, and\n", "minimising that at fixed $\\int t\\,\\mathrm{d}x$ gives $t \\propto \\sqrt{M}$; under uniform\n", "load $M \\propto (L-x)^2$, so the optimum is the linear wedge $p = 1$. That is a real\n", "answer to recover, so we start $p$ deliberately away from it." ] }, { "cell_type": "code", "execution_count": 2, "id": "7113c623", "metadata": { "execution": { "iopub.execute_input": "2026-09-09T21:40:21.752081Z", "iopub.status.busy": "2026-09-09T21:40:21.751892Z", "iopub.status.idle": "2026-09-09T21:40:21.943981Z", "shell.execute_reply": "2026-09-09T21:40:21.943641Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "starting point: p = 0.3, compliance 1.689510e-01, mass 3.0465\n" ] } ], "source": [ "recorder = csdl.Recorder(inline=True)\n", "recorder.start()\n", "\n", "domain = hm.ShellDomain(mesh, element=\"CG2CG1\")\n", "bcs = hm.clamp(domain, where=hm.near(\"x\", 0.0))\n", "loads = hm.pressure(domain, PRESSURE_Z)\n", "\n", "EPS = 1e-3 # keeps the tip thickness positive\n", "xi = (LENGTH - domain.node_coords[:, 0]) / LENGTH + EPS\n", "\n", "# t(x) = A * xi**p. A variable exponent needs exp(p log xi); xi is a constant array,\n", "# so this stays a clean CSDL expression in the two design variables.\n", "scale = csdl.Variable(value=np.array([H_VAL]), name=\"scale\")\n", "exponent = csdl.Variable(value=np.array([0.3]), name=\"exponent\") # away from p = 1\n", "log_xi = csdl.Variable(value=np.log(xi))\n", "thickness = csdl.expand(scale, xi.shape) * csdl.exp(csdl.expand(exponent, xi.shape) * log_xi)\n", "\n", "material = hm.isotropic(domain, E=E_VAL, nu=NU_VAL,\n", " thickness=hm.from_nodal(domain, thickness), density=RHO_VAL)\n", "state = hm.solve(domain, material, loads, bcs)\n", "compliance, mass = hm.compliance(state), hm.mass(state)\n", "\n", "c0 = float(np.ravel(compliance.value)[0])\n", "mass_target = RHO_VAL * H_VAL * LENGTH * WIDTH\n", "print(f\"starting point: p = 0.3, compliance {c0:.6e}, mass {float(np.ravel(mass.value)[0]):.4f}\")" ] }, { "cell_type": "code", "execution_count": 3, "id": "511453be", "metadata": { "execution": { "iopub.execute_input": "2026-09-09T21:40:21.945221Z", "iopub.status.busy": "2026-09-09T21:40:21.945129Z", "iopub.status.idle": "2026-09-09T21:40:23.871116Z", "shell.execute_reply": "2026-09-09T21:40:23.870771Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Optimization terminated successfully (Exit mode 0)\n", " Final objective value : 4.205578e-02\n", " Final optimality : 1.892155e-14\n", " Final feasibility : 3.002043e-13\n", " Number of major iterations : 8\n", " Number of function evaluations : 8\n", " Number of derivative evaluations : 8\n", " Average Function evaluation time : 0.092107 s per evaluation\n", " Average Derivative evaluation time : 0.077418 s per evaluation\n", " Total Function evaluation time : 0.736854 s [ 54.26%]\n", " Total Derivative evaluation time : 0.619346 s [ 45.61%]\n", " Optimizer time : 0.000200 s [ 0.01%]\n", " Processing time : 0.001602 s [ 0.12%]\n", " Visualization time : 0.000000 s [ 0.00%]\n", " Total optimization time : 1.358003 s [100.00%]\n", " Summary saved to : tutorial_taper_outputs/2026-09-09_14.40.22.510003/slsqp_summary.out\n" ] } ], "source": [ "scale.set_as_design_variable(lower=1e-3, upper=1.0)\n", "exponent.set_as_design_variable(lower=0.0, upper=3.0)\n", "mass.set_as_constraint(lower=mass_target, upper=mass_target)\n", "compliance.set_as_objective()\n", "\n", "from modopt import CSDLAlphaProblem, PySLSQP\n", "\n", "sim = csdl.experimental.PySimulator(recorder)\n", "problem = CSDLAlphaProblem(problem_name=\"tutorial_taper\", simulator=sim)\n", "PySLSQP(problem, solver_options={\"maxiter\": 200, \"acc\": 1e-12}).solve()\n", "\n", "recorder.stop()" ] }, { "cell_type": "code", "execution_count": 4, "id": "8efa7ba2", "metadata": { "execution": { "iopub.execute_input": "2026-09-09T21:40:23.872325Z", "iopub.status.busy": "2026-09-09T21:40:23.872090Z", "iopub.status.idle": "2026-09-09T21:40:23.874556Z", "shell.execute_reply": "2026-09-09T21:40:23.874217Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "recovered exponent : p = 1.0603 (analytic 1.0)\n", "compliance : 1.689510e-01 -> 4.205578e-02 (75.1% lower)\n", "mass : 4.000000 (target 4.0)\n" ] } ], "source": [ "p_opt = float(np.ravel(exponent.value)[0])\n", "c1 = float(np.ravel(compliance.value)[0])\n", "\n", "print(f\"recovered exponent : p = {p_opt:.4f} (analytic 1.0)\")\n", "print(f\"compliance : {c0:.6e} -> {c1:.6e} ({100 * (1 - c1 / c0):.1f}% lower)\")\n", "print(f\"mass : {float(np.ravel(mass.value)[0]):.6f} (target {mass_target})\")" ] }, { "cell_type": "markdown", "id": "8d37a22f", "metadata": {}, "source": [ "The optimizer recovers the analytic exponent to within about 6% on this deliberately\n", "tiny mesh, driving through the shell solve's adjoint from a starting point well away\n", "from the answer. The remaining error is honest: the plate has finite width and\n", "transverse shear, and the tip thickness is $\\varepsilon$ rather than zero, so it is not\n", "exactly the beam the closed form assumes. `ex_optimal_thickness_taper.py` runs this as\n", "a gated benchmark on finer meshes and holds it to 5%.\n", "\n", "That benchmark's docstring also records the measurements behind the ill-posedness\n", "warning above \u2014 on its fixture, at equal mass, uniform gives $C = 0.056$, the analytic\n", "wedge $C = 0.017$, and a free-nodal optimum $C = 0.0037$. The checkerboard beats the\n", "wedge by five times, and it is *right* to.\n", "`examples/advanced_examples/ex_thickness_opt.py` runs the free-nodal problem if you\n", "want to see that for yourself." ] }, { "cell_type": "markdown", "id": "fa70e42f", "metadata": {}, "source": [ "## Shape (mesh-coordinate) derivatives\n", "\n", "Hermit differentiates the residual and every output functional with respect to\n", "`ufl.SpatialCoordinate`, then scatters the sensitivity back to your node ordering.\n", "This is the gradient a shape optimizer, or a chain rule back to a geometry\n", "parameterisation, consumes.\n", "\n", "Make the geometry a live input by passing `geometry=` a `hm.geometry(...)` built from a\n", "`csdl.Variable` \u2014 either `node_disp=` (a perturbation added to the reference\n", "coordinates) or `nodes=` (absolute coordinates). Leave it out and the geometry is a\n", "constant: no mesh-derivative forms are built, so there is no cost to not using this." ] }, { "cell_type": "code", "execution_count": 5, "id": "5ac6d8e4", "metadata": { "execution": { "iopub.execute_input": "2026-09-09T21:40:23.875387Z", "iopub.status.busy": "2026-09-09T21:40:23.875271Z", "iopub.status.idle": "2026-09-09T21:40:24.156078Z", "shell.execute_reply": "2026-09-09T21:40:24.155728Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "shape-derivative field: (33, 3) (one gradient vector per mesh node)\n" ] } ], "source": [ "recorder = csdl.Recorder(inline=True)\n", "recorder.start()\n", "\n", "domain = hm.ShellDomain(plate_mesh(nx, ny), element=\"CG2CG1\")\n", "bcs = hm.clamp(domain, where=hm.near(\"x\", 0.0))\n", "material = hm.isotropic(domain, E=E_VAL, nu=NU_VAL, thickness=H_VAL, density=RHO_VAL)\n", "\n", "node_disp = csdl.Variable(value=np.zeros((domain.n_nodes, 3)), name=\"node_disp\")\n", "state = hm.solve(domain, material, hm.pressure(domain, PRESSURE_Z), bcs,\n", " geometry=hm.geometry(domain, node_disp=node_disp))\n", "compliance = hm.compliance(state)\n", "\n", "sim = csdl.experimental.PySimulator(recorder)\n", "dC_dX = np.asarray(\n", " sim.compute_totals([compliance], [node_disp])[compliance, node_disp]\n", ").reshape(domain.n_nodes, 3)\n", "\n", "recorder.stop()\n", "\n", "print(f\"shape-derivative field: {dC_dX.shape} (one gradient vector per mesh node)\")" ] }, { "cell_type": "markdown", "id": "ef09269f", "metadata": {}, "source": [ "`dC_dX[k]` is $\\partial(\\text{compliance})/\\partial\\mathbf{x}_k$ at mesh node $k$, in\n", "your node ordering. Two checks worth doing on a new model:\n", "\n", "Stretching the plate along its axis makes it more compliant, so $\\partial C/\\partial x$\n", "should be positive at the tip and negative at the root \u2014 moving material away from the\n", "support in one case, towards it in the other.\n", "\n", "The out-of-plane channel should be **exactly** zero on a flat plate: perturbing a node\n", "by $+\\delta z$ and by $-\\delta z$ gives mirror-image configurations with identical\n", "compliance, so the first-order term vanishes by symmetry. That is a genuine result, not\n", "a dead channel \u2014 give the plate a shallow curve and the same component is nonzero." ] }, { "cell_type": "code", "execution_count": 6, "id": "e6cfa71c", "metadata": { "execution": { "iopub.execute_input": "2026-09-09T21:40:24.156915Z", "iopub.status.busy": "2026-09-09T21:40:24.156812Z", "iopub.status.idle": "2026-09-09T21:40:24.159181Z", "shell.execute_reply": "2026-09-09T21:40:24.158857Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "mean dC/dx at the tip : +2.3078e-02\n", "mean dC/dx at the root : -2.2939e-02\n", "max |dC/dz|, flat plate: 0.0000e+00 (zero by reflection symmetry)\n" ] } ], "source": [ "x = domain.node_coords[:, 0]\n", "tip, root = x > LENGTH - 1e-9, x < 1e-9\n", "print(f\"mean dC/dx at the tip : {dC_dX[tip, 0].mean():+.4e}\")\n", "print(f\"mean dC/dx at the root : {dC_dX[root, 0].mean():+.4e}\")\n", "print(f\"max |dC/dz|, flat plate: {np.abs(dC_dX[:, 2]).max():.4e} (zero by reflection symmetry)\")" ] }, { "cell_type": "code", "execution_count": 7, "id": "85ff5595", "metadata": { "execution": { "iopub.execute_input": "2026-09-09T21:40:24.159898Z", "iopub.status.busy": "2026-09-09T21:40:24.159804Z", "iopub.status.idle": "2026-09-09T21:40:24.444693Z", "shell.execute_reply": "2026-09-09T21:40:24.444329Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "max |dC/dz|, cambered : 2.0848e-03\n" ] } ], "source": [ "# the same gradient about a shallow parabolic camber -- the z channel is live\n", "recorder = csdl.Recorder(inline=True)\n", "recorder.start()\n", "\n", "curved = hm.ShellDomain(plate_mesh(nx, ny), element=\"CG2CG1\")\n", "camber = np.zeros((curved.n_nodes, 3))\n", "camber[:, 2] = 0.5 * (curved.node_coords[:, 0] / LENGTH) ** 2\n", "nd = csdl.Variable(value=camber, name=\"node_disp\")\n", "\n", "curved_state = hm.solve(\n", " curved,\n", " hm.isotropic(curved, E=E_VAL, nu=NU_VAL, thickness=H_VAL, density=RHO_VAL),\n", " hm.pressure(curved, PRESSURE_Z),\n", " hm.clamp(curved, where=hm.near(\"x\", 0.0)),\n", " geometry=hm.geometry(curved, node_disp=nd),\n", ")\n", "C_curved = hm.compliance(curved_state)\n", "\n", "sim = csdl.experimental.PySimulator(recorder)\n", "g = np.asarray(sim.compute_totals([C_curved], [nd])[C_curved, nd]).reshape(curved.n_nodes, 3)\n", "recorder.stop()\n", "\n", "print(f\"max |dC/dz|, cambered : {np.abs(g[:, 2]).max():.4e}\")" ] }, { "cell_type": "markdown", "id": "ca80cb1c", "metadata": {}, "source": [ "### Requirements and caveats\n", "\n", "- **Serial mesh** whose `input_global_indices` form a node permutation. `ShellDomain`\n", " enforces this for every use, not only for shape derivatives.\n", "- **UFL patch.** UFL's `CoordinateDerivative` handler crashes on the shell bending\n", " term, in both UFL 2024.2 (DOLFINx 0.9) and 2026.1 (DOLFINx 0.11). Hermit applies a\n", " small runtime monkeypatch on `import hermit`; the installed UFL files are untouched.\n", "- **Non-smooth stabilization scale.** UFL drops the non-differentiable subgradient of\n", " `CellDiameter`, which appears only in the drilling / penalty stabilization scale. The\n", " effect on outputs is at the 1e-5 level.\n", "- **Penalty conditioning.** Under a very large geometry perturbation the penalty-BC\n", " system becomes ill-conditioned and finite-difference checks stop converging. That is\n", " a conditioning limit, not a bug \u2014 use the strong-BC path or a smaller step to verify.\n", "\n", "The mesh-coordinate adjoints are finite-difference validated in\n", "`tests/test_mesh_coord_deriv.py`, which measures a few parts in $10^{5}$ relative\n", "error and gates at $3\\times10^{-3}$; `ex_shape_derivative_fd.py` gates at\n", "$5\\times10^{-3}$. Those gates are set by finite-difference truncation and\n", "cancellation, not by the adjoint." ] } ], "metadata": { "kernelspec": { "display_name": "hermit", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.13" } }, "nbformat": 4, "nbformat_minor": 5 }