|
70 | 70 | "metadata": {}, |
71 | 71 | "outputs": [], |
72 | 72 | "source": [ |
73 | | - "# The PyMAPDL Host (X.X.X.X IP Address) and Port (XXX) are pulled from the environment \n", |
| 73 | + "# The PyMAPDL Host (X.X.X.X IP Address) and Port (XXX) are pulled from the environment\n", |
74 | 74 | "import os\n", |
| 75 | + "\n", |
75 | 76 | "HOST = os.getenv(\"MAPDL_HOST\")\n", |
76 | 77 | "if HOST is None:\n", |
77 | 78 | " raise ValueError(\n", |
|
83 | 84 | " raise ValueError(\n", |
84 | 85 | " \"Unable to read $MAPDL_PORT from the environment. \"\n", |
85 | 86 | " + \"Use 'export MAPDL_PORT=X' for the port of your MAPDL Instance.\"\n", |
86 | | - " )\n" |
| 87 | + " )" |
87 | 88 | ] |
88 | 89 | }, |
89 | 90 | { |
|
112 | 113 | "import matplotlib.pyplot as plt\n", |
113 | 114 | "from tesseract_core import Tesseract\n", |
114 | 115 | "from tesseract_core.runtime.experimental import TesseractReference\n", |
115 | | - "from utils import hex_to_pyvista, plot_grid, plot_grid_slice, plot_mesh, stop_grads_int\n", |
| 116 | + "from utils import hex_to_pyvista, plot_grid, plot_grid_slice, plot_mesh\n", |
116 | 117 | "\n", |
117 | 118 | "from tesseract_jax import apply_tesseract" |
118 | 119 | ] |
|
466 | 467 | } |
467 | 468 | ], |
468 | 469 | "source": [ |
469 | | - "import pyvista as pv\n", |
470 | 470 | "import numpy as np\n", |
| 471 | + "import pyvista as pv\n", |
471 | 472 | "\n", |
472 | 473 | "grid = pv.ImageData(\n", |
473 | 474 | " dimensions=np.array((Nx, Ny, Nz)) + 1,\n", |
474 | | - " origin=(-Lx / 2, -Ly / 2, -Lz / 2),\n", |
| 475 | + " origin=(-Lx / 2, -Ly / 2, -Lz / 2),\n", |
475 | 476 | " spacing=(Lx / Nx, Ly / Ny, Lz / Nz), # These are the cell sizes along each axis\n", |
476 | 477 | ")\n", |
477 | 478 | "# repeated casts will eventually expose cell_connectivitiy\n", |
|
541 | 542 | }, |
542 | 543 | { |
543 | 544 | "cell_type": "code", |
544 | | - "execution_count": 13, |
| 545 | + "execution_count": null, |
545 | 546 | "id": "9ea51a40", |
546 | 547 | "metadata": { |
547 | 548 | "editable": true, |
|
572 | 573 | " Args:\n", |
573 | 574 | " mesh: Dictionary containing 'points' array.\n", |
574 | 575 | " Lx: Size of the domain in x-direction.\n", |
| 576 | + " Ly: Size of the domain in y-direction.\n", |
575 | 577 | " Lz: Size of the domain in z-direction.\n", |
576 | 578 | "\n", |
577 | 579 | " Returns:\n", |
|
582 | 584 | "\n", |
583 | 585 | " dirichlet_mask = pts[:, 0] <= -Lx / 2\n", |
584 | 586 | " # von Neumann condition (select nodes at x=Lx with constraints on y and z)\n", |
585 | | - " x_lim = Lx/2\n", |
| 587 | + " x_lim = Lx / 2\n", |
586 | 588 | " y_min = -0.1 * Ly\n", |
587 | 589 | " y_max = 0.1 * Ly\n", |
588 | 590 | " z_min = -0.1 * Lz\n", |
|
599 | 601 | " ),\n", |
600 | 602 | " )\n", |
601 | 603 | "\n", |
602 | | - "\n", |
603 | 604 | " return dirichlet_mask, van_neumann_mask\n", |
604 | 605 | "\n", |
605 | 606 | "\n", |
|
653 | 654 | "outputs": [], |
654 | 655 | "source": [ |
655 | 656 | "# load the Finite Element Method Tesseract\n", |
656 | | - "fem_tess = Tesseract.from_tesseract_api(\"pymapdl_tess/tesseract_api.py\")\n" |
| 657 | + "fem_tess = Tesseract.from_tesseract_api(\"pymapdl_tess/tesseract_api.py\")" |
657 | 658 | ] |
658 | 659 | }, |
659 | 660 | { |
|
877 | 878 | }, |
878 | 879 | { |
879 | 880 | "cell_type": "code", |
880 | | - "execution_count": 17, |
| 881 | + "execution_count": null, |
881 | 882 | "id": "6c1dae2f", |
882 | 883 | "metadata": {}, |
883 | 884 | "outputs": [], |
|
942 | 943 | "\n", |
943 | 944 | " van_neumann_mask = jax.lax.stop_gradient(van_neumann_mask)\n", |
944 | 945 | " dirichlet_mask = jax.lax.stop_gradient(dirichlet_mask)\n", |
945 | | - " dirichlet_values = jnp.array([0.0])\n", |
946 | | - " van_neumann_values = jnp.array([[0.0, 0.0, 0.1]])\n", |
| 946 | + " # dirichlet_values = jnp.array([0.0])\n", |
| 947 | + " # van_neumann_values = jnp.array([[0.0, 0.0, 0.1]])\n", |
947 | 948 | "\n", |
948 | | - " # Instead of passing all inputs and trying to stop_gradient on them,\n", |
949 | | - " # we need to wrap the tesseract call to only allow gradients w.r.t. rho\n", |
950 | 949 | " # TODO unify this framework\n", |
951 | 950 | " dirichlet_idx = np.where(dirichlet_mask)[0]\n", |
952 | 951 | " von_neumann_idx = np.where(van_neumann_mask)[0]\n", |
953 | | - " von_neumann_values = np.array([0, 0.0, 0.1]) + np.zeros((von_neumann_idx.shape[0], 3))\n", |
954 | | - " \n", |
| 952 | + " von_neumann_values = np.array([0, 0.0, 0.1]) + np.zeros(\n", |
| 953 | + " (von_neumann_idx.shape[0], 3)\n", |
| 954 | + " )\n", |
| 955 | + "\n", |
955 | 956 | " # TODO again note the transpose on rho!\n", |
956 | 957 | " rho_for_fea = rho_grid.T.reshape((hex_mesh[\"faces\"].shape[0], 1))\n", |
957 | 958 | "\n", |
|
973 | 974 | " },\n", |
974 | 975 | " )[\"compliance\"]\n", |
975 | 976 | "\n", |
976 | | - " return c, {\n", |
977 | | - " \"sdf\": sdf,\n", |
978 | | - " \"rho_grid\": rho_grid\n", |
979 | | - " }\n", |
| 977 | + " return c, {\"sdf\": sdf, \"rho_grid\": rho_grid}\n", |
980 | 978 | "\n", |
981 | 979 | "\n", |
982 | 980 | "identity_and_store_grads.defvjp(identity_fwd, identity_bwd)\n", |
983 | | - "grad_fn = jax.value_and_grad(loss, has_aux=True)\n" |
| 981 | + "grad_fn = jax.value_and_grad(loss, has_aux=True)" |
984 | 982 | ] |
985 | 983 | }, |
986 | 984 | { |
|
1200 | 1198 | "for i in range(num_tests):\n", |
1201 | 1199 | " print(i)\n", |
1202 | 1200 | " params_up = params.at[i].add(FD_delta)\n", |
1203 | | - " (fupp, _), _ = grad_fn(params_up, iteration=i*2)\n", |
| 1201 | + " (fupp, _), _ = grad_fn(params_up, iteration=i * 2)\n", |
1204 | 1202 | "\n", |
1205 | 1203 | " if run_central_difference:\n", |
1206 | 1204 | " params_down = params.at[i].subtract(2.0 * FD_delta)\n", |
1207 | | - " (fdown, _), _ = grad_fn(params_down, iteration=i*2+1)\n", |
| 1205 | + " (fdown, _), _ = grad_fn(params_down, iteration=i * 2 + 1)\n", |
1208 | 1206 | " FD_grads[i] = (fupp - fdown) / FD_delta / 2.0\n", |
1209 | 1207 | "\n", |
1210 | 1208 | " else:\n", |
|
1235 | 1233 | }, |
1236 | 1234 | { |
1237 | 1235 | "cell_type": "code", |
1238 | | - "execution_count": 18, |
| 1236 | + "execution_count": null, |
1239 | 1237 | "id": "849a7559-90fc-4531-a0ed-54bf4a1c7817", |
1240 | 1238 | "metadata": { |
1241 | 1239 | "editable": true, |
|
1248 | 1246 | "source": [ |
1249 | 1247 | "import mmapy\n", |
1250 | 1248 | "import numpy as np\n", |
1251 | | - "from tesseract_core.runtime import Array, Differentiable, Float32, Int32\n", |
| 1249 | + "from tesseract_core.runtime import Array, Float32, Int32\n", |
1252 | 1250 | "\n", |
1253 | 1251 | "\n", |
1254 | 1252 | "class MMAOptimizer:\n", |
1255 | 1253 | " \"\"\"A wrapper for the MMA optimizer from mmapy.\n", |
1256 | | - " github.com/arjendeetman/GCMMA-MMA-Python\n", |
| 1254 | + "\n", |
| 1255 | + " Source is github.com/arjendeetman/GCMMA-MMA-Python.\n", |
1257 | 1256 | " mmapy is a pretty barebones implementation of MMA in python. It should work for now.\n", |
1258 | 1257 | " Alternatives to consider:\n", |
1259 | 1258 | " - github.com/LLNL/pyMMAopt\n", |
1260 | | - " - pyopt.org/reference/optimizers.mma.html\n", |
| 1259 | + " - pyopt.org/reference/optimizers.mma.html.\n", |
1261 | 1260 | "\n", |
1262 | 1261 | " \"\"\"\n", |
1263 | 1262 | "\n", |
|
1317 | 1316 | " if iteration == 1:\n", |
1318 | 1317 | " self.objective_scale_factor = np.abs(self.objective_scale / objective_value)\n", |
1319 | 1318 | " objective_value *= self.objective_scale_factor\n", |
1320 | | - " objective_gradient = np.asarray(objective_gradient) * self.objective_scale_factor\n", |
| 1319 | + " objective_gradient = (\n", |
| 1320 | + " np.asarray(objective_gradient) * self.objective_scale_factor\n", |
| 1321 | + " )\n", |
1321 | 1322 | "\n", |
1322 | 1323 | " # the bounds dont necessarily change every iteration\n", |
1323 | 1324 | " if x_min is None:\n", |
|
1335 | 1336 | " )\n", |
1336 | 1337 | "\n", |
1337 | 1338 | " # calculate the next iteration of x\n", |
1338 | | - " xmma, ymma, zmma, lam, xsi, eta, mu, zet, s, low, upp = mmapy.mmasub(\n", |
| 1339 | + " xmma, _ymma, _zmma, _lam, _xsi, _eta, _mu, _zet, _s, low, upp = mmapy.mmasub(\n", |
1339 | 1340 | " self.m,\n", |
1340 | 1341 | " self.n,\n", |
1341 | 1342 | " iteration,\n", |
|
1378 | 1379 | " constraint_gradients: Array[(None, None), Float32] = None,\n", |
1379 | 1380 | " ) -> None:\n", |
1380 | 1381 | " def check_shape(shape: tuple, expected_shape: tuple, name: str) -> None:\n", |
1381 | | - " if (len(shape)==1) or (shape[0] != expected_shape[0] or shape[1] != expected_shape[1]):\n", |
| 1382 | + " if (len(shape) == 1) or (\n", |
| 1383 | + " shape[0] != expected_shape[0] or shape[1] != expected_shape[1]\n", |
| 1384 | + " ):\n", |
1382 | 1385 | " raise TypeError(\n", |
1383 | 1386 | " f\"MMAError: The '{name}' was expected to have shape {expected_shape} but has shape {shape}.\"\n", |
1384 | 1387 | " )\n", |
|
1393 | 1396 | " if constraint_gradients is not None:\n", |
1394 | 1397 | " check_shape(\n", |
1395 | 1398 | " constraint_gradients.shape, (self.m, self.n), \"constraint gradients\"\n", |
1396 | | - " )\n" |
| 1399 | + " )" |
1397 | 1400 | ] |
1398 | 1401 | }, |
1399 | 1402 | { |
1400 | 1403 | "cell_type": "code", |
1401 | | - "execution_count": 19, |
| 1404 | + "execution_count": null, |
1402 | 1405 | "id": "39e30381", |
1403 | 1406 | "metadata": { |
1404 | 1407 | "editable": true, |
|
1739 | 1742 | "delta_x = 0.0\n", |
1740 | 1743 | "delta_y = 0.0\n", |
1741 | 1744 | "delta_z = Lz / 3\n", |
1742 | | - "eps = 1.e-3 # a small value to ease numerics\n", |
| 1745 | + "eps = 1.0e-3 # a small value to ease numerics\n", |
1743 | 1746 | "bar_params = np.array(initial_params).reshape(n_chains, n_edges_per_chain + 1, 3)\n", |
1744 | 1747 | "param_min = bar_params - eps\n", |
1745 | 1748 | "param_max = bar_params + eps\n", |
1746 | 1749 | "for bar_idx, bar in enumerate(bar_params):\n", |
1747 | | - " for coord_idx, xyz in enumerate(bar):\n", |
| 1750 | + " for coord_idx, _xyz in enumerate(bar):\n", |
1748 | 1751 | " for i, delta in enumerate([delta_x, delta_y, delta_z]):\n", |
1749 | 1752 | " param_min[bar_idx][coord_idx][i] -= delta\n", |
1750 | 1753 | " param_max[bar_idx][coord_idx][i] += delta\n", |
|
1762 | 1765 | "loss_hist = []\n", |
1763 | 1766 | "params_hist = []\n", |
1764 | 1767 | "aux_hist = []\n", |
1765 | | - "optimizer = MMAOptimizer(initial_params[:,None], param_min, param_max, num_constraints, x_update_limit=x_update_limit)\n", |
| 1768 | + "optimizer = MMAOptimizer(\n", |
| 1769 | + " initial_params[:, None],\n", |
| 1770 | + " param_min,\n", |
| 1771 | + " param_max,\n", |
| 1772 | + " num_constraints,\n", |
| 1773 | + " x_update_limit=x_update_limit,\n", |
| 1774 | + ")\n", |
1766 | 1775 | "\n", |
1767 | 1776 | "\n", |
1768 | 1777 | "for i in range(n_steps):\n", |
1769 | 1778 | " (loss_value, aux), grads = grad_fn(params, iteration=i)\n", |
1770 | 1779 | " print(loss_value)\n", |
1771 | 1780 | " print(params)\n", |
1772 | 1781 | " print(grads)\n", |
1773 | | - " np_grads = np.array(grads[:,None])\n", |
1774 | | - " np_params = np.array(params[:,None])\n", |
1775 | | - " np_params = optimizer.calculate_next_x(loss_value, np_grads, g, dgdx, i+1, np_params)\n", |
| 1782 | + " np_grads = np.array(grads[:, None])\n", |
| 1783 | + " np_params = np.array(params[:, None])\n", |
| 1784 | + " np_params = optimizer.calculate_next_x(\n", |
| 1785 | + " loss_value, np_grads, g, dgdx, i + 1, np_params\n", |
| 1786 | + " )\n", |
1776 | 1787 | " params = jnp.array(np_params.flatten())\n", |
1777 | | - " \n", |
| 1788 | + "\n", |
1778 | 1789 | " loss_hist.append(loss_value)\n", |
1779 | 1790 | " params_hist.append(params)\n", |
1780 | 1791 | " aux_hist.append(aux)\n", |
1781 | 1792 | "\n", |
1782 | | - "\n", |
1783 | 1793 | " print(f\"Iteration {i + 1}, Loss: {loss_value:.4f}\")" |
1784 | 1794 | ] |
1785 | 1795 | }, |
|
1845 | 1855 | " mesh = hex_mesh\n", |
1846 | 1856 | " rho_dot = grad_storage[i + 2000][1][: len(mesh[\"faces\"])]\n", |
1847 | 1857 | " pv_mesh = hex_to_pyvista(\n", |
1848 | | - " mesh[\"points\"], mesh[\"faces\"], {\"rho\": aux_hist[i][\"rho_grid\"].T.flatten(), \"rho_dot\": rho_dot.T.flatten()}\n", |
| 1858 | + " mesh[\"points\"],\n", |
| 1859 | + " mesh[\"faces\"],\n", |
| 1860 | + " {\"rho\": aux_hist[i][\"rho_grid\"].T.flatten(), \"rho_dot\": rho_dot.T.flatten()},\n", |
1849 | 1861 | " )\n", |
1850 | 1862 | " pv_mesh.save(f\"vtks/fem_shapeopt_mesh{i + 1}.vtk\")" |
1851 | 1863 | ] |
|
1882 | 1894 | ], |
1883 | 1895 | "metadata": { |
1884 | 1896 | "kernelspec": { |
1885 | | - "display_name": "tesseract-jax", |
| 1897 | + "display_name": "fem", |
1886 | 1898 | "language": "python", |
1887 | | - "name": "tesseract-jax" |
| 1899 | + "name": "python3" |
1888 | 1900 | }, |
1889 | 1901 | "language_info": { |
1890 | 1902 | "codemirror_mode": { |
|
1896 | 1908 | "name": "python", |
1897 | 1909 | "nbconvert_exporter": "python", |
1898 | 1910 | "pygments_lexer": "ipython3", |
1899 | | - "version": "3.13.5" |
| 1911 | + "version": "3.13.9" |
1900 | 1912 | } |
1901 | 1913 | }, |
1902 | 1914 | "nbformat": 4, |
|
0 commit comments