Skip to content

Commit

Permalink
Make Toni's changes
Browse files Browse the repository at this point in the history
  • Loading branch information
@donaldcummins
donaldcummins committed Dec 19, 2024
1 parent 565205c commit a4874e8
Showing 1 changed file with 46 additions and 37 deletions.
83 changes: 46 additions & 37 deletions Design_Optimisation_Metamodels.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -88,8 +88,8 @@
"\n",
"$$\n",
"\\begin{aligned}\n",
"u = u_{in} &= 6 \\omega (1-x_2)x_2 \\boldsymbol{e}_1 \\quad &\\textrm{ on } \\Gamma_{in} \\\\\n",
"u = u_{bc} &= \\frac{2(1-\\omega)}{\\tan \\theta} \\boldsymbol{e}_1 - 2(1-\\omega) \\boldsymbol{e}_2 \\quad &\\textrm{ on } \\Gamma_{bc}\\\\\n",
"u &= 6 \\: v_{max} \\: \\omega (1-x_2)x_2 \\boldsymbol{e}_1 \\quad &\\textrm{ on } \\Gamma_{in} \\\\\n",
"u &= 6 \\: v_{max} \\: (1-\\omega)(2-x_1)(x_1-1)[ \\tan^{-1}(\\theta) \\boldsymbol{e}_1 - \\boldsymbol{e}_2] \\quad &\\textrm{ on } \\Gamma_{bc}\\\\\n",
"p &= 0 &\\textrm{ on } \\Gamma_{out} \\\\\n",
"u &= 0 &\\textrm{ on } \\Gamma_{w}.\n",
"\\end{aligned}\n",
Expand Down Expand Up @@ -279,7 +279,7 @@
" # Define physical parameters of flow problem\n",
" mu = 0.035 # dynamic viscosity (value for blood)\n",
" rho = 1.060 # density (value for blood)\n",
" Re = 90 # desired Reynolds number for pipe flow (biological Re is between 100-200 in coronary arteries)\n",
" Re = 100 # desired Reynolds number for pipe flow (biological Re is between 100-200 in coronary arteries)\n",
"\n",
" # Define computational domain\n",
" L = 5.0\n",
Expand All @@ -299,10 +299,10 @@
" # Define the subdomain \\Omega_obs and its complement\n",
" class OmegaObs(SubDomain):\n",
" def inside(self, x, on_boundary):\n",
" return True if (x[0] >= 1 and x[0] <= 4) else False\n",
" return True if (x[0] >= 1.0 and x[0] <= 4.0) else False\n",
" class OmegaRest(SubDomain):\n",
" def inside(self, x, on_boundary):\n",
" return True if (x[0] <= 1 or x[0] >= 4) else False\n",
" return True if (x[0] <= 1.0 or x[0] >= 4.0) else False\n",
"\n",
" subdomains = MeshFunction('size_t', mesh, mesh.topology().dim())\n",
" subdomain0 = OmegaObs()\n",
Expand All @@ -323,17 +323,17 @@
" inflow = 'near(x[0], 0)'\n",
" outflow = 'near(x[0], 5.0)'\n",
" bottom = 'near(x[1], 0)'\n",
" top = 'near(x[1], 1) && ((x[0]<1.0) || (x[0]>1.5))'\n",
" anastomosis = 'near(x[1], 1) && ((x[0]>1.0) && (x[0]<1.5))'\n",
" top = 'near(x[1], 1) && ((x[0]<1.0) || (x[0]>2.0))'\n",
" anastomosis = 'near(x[1], 1) && ((x[0]>1.0) && (x[0]<2.0))'\n",
"\n",
" # Scale flow appropriately to the Reynolds number\n",
" ReMultiplier = Re / (rho/mu)\n",
" VMax = Re / (rho/mu)\n",
"\n",
" # Define parabolic inflow profile\n",
" inflow_profile = (str(omega * ReMultiplier) + '*6*x[1]*(1-x[1])', '0')\n",
" inflow_profile = (str(omega * VMax) + '*6*x[1]*(1-x[1])', '0')\n",
"\n",
" # Define anastomosis inflow profile parameterised by angle\n",
" anastomosis_profile = (str((1-omega)*ReMultiplier) +'*2/tan(' + str(theta * 2*pi/360) + ')', '-2*' + str((1-omega)*ReMultiplier))\n",
" anastomosis_profile = (str((1-omega)*VMax) + ' * (x[0]-1.0) * (2.0-x[0]) * 6 /tan(' + str(theta*100 * 2*pi/360) + ')', str((1-omega)*VMax) + ' * (-1.0) * (x[0]-1.0) * (2.0-x[0]) * 6' )\n",
"\n",
" # Define boundary conditions\n",
" bcu_inflow = DirichletBC(W.sub(0), Expression(inflow_profile, degree=2), inflow)\n",
Expand Down Expand Up @@ -409,7 +409,7 @@
" # Solve Navier-Stokes and evaluate cost functional\n",
" E, V = solve_navier_stokes(W, nu, bcs, IF_OmegaObs)\n",
"\n",
" return V"
" return E"
]
},
{
Expand All @@ -432,7 +432,7 @@
},
"outputs": [],
"source": [
"solveAnastomosisProblemNL(60, 0.3)"
"solveAnastomosisProblemNL(0.4, 0.3)"
]
},
{
Expand Down Expand Up @@ -543,7 +543,7 @@
},
"outputs": [],
"source": [
"Ns = 49 # Number of snapshots to generate\n",
"Ns = 113 # Number of snapshots to generate\n",
"\n",
"snapshots = np.empty((Ns, 3))\n",
"\n",
Expand All @@ -555,9 +555,9 @@
"\n",
" case \"SmolyakInterpolation\":\n",
" dim = 2\n",
" level = 3\n",
" level = 4\n",
" # Use Clenshaw-Curtis points including the boundary on [0,1]^2\n",
" grid = pysgpp.Grid.createLinearClenshawCurtisBoundaryGrid(dim)\n",
" grid = pysgpp.Grid.createPolyClenshawCurtisBoundaryGrid(dim,degree=2)\n",
" gridStorage = grid.getStorage()\n",
" grid.getGenerator().regular(level)\n",
" print(\"number of grid points: {}\".format(gridStorage.getSize()))\n",
Expand All @@ -568,20 +568,20 @@
" case \"RBFInterpolation/RegularGrid\":\n",
" idx = float(n % 7) / 6.0\n",
" jdx = float(floor(float(n / 7.0))) / 6.0\n",
" theta = idx * 40. + 30.\n",
" theta = idx * 0.4 + 0.3\n",
" omega = jdx * 0.5 + 0.3\n",
"\n",
" case \"RBFInterpolation/RandomPoints\":\n",
" theta = random.random() * 40 + 30 # Create random angle between 30 and 70 degrees\n",
" theta = random.random() * 0.4 + 30 # Create random angle between 30 and 70 degrees\n",
" omega = random.random() * 0.5 + 0.3 # Create random flow split between 30% and 80%\n",
"\n",
" case \"RBFInterpolation/HaltonSequence\":\n",
" theta = HaltonSample[n,0] * 40 + 30\n",
" theta = HaltonSample[n,0] * 0.4 + 0.3\n",
" omega = HaltonSample[n,1] * 0.5 + 0.3\n",
"\n",
" case \"SmolyakInterpolation\":\n",
" gp = gridStorage.getPoint(n)\n",
" theta = gp.getStandardCoordinate(0) * 40 + 30\n",
" theta = gp.getStandardCoordinate(0) * 0.4 + 0.3\n",
" omega = gp.getStandardCoordinate(1) * 0.5 + 0.3\n",
"\n",
" # Solve the PDE for the values (theta,omega)\n",
Expand Down Expand Up @@ -614,7 +614,10 @@
" if (x.ndim != 2):\n",
" x = x.reshape(1,2)\n",
"\n",
" return (RBFInterpolator(snapshots[:,0:2:1], snapshots[:,2], kernel='inverse_multiquadric', epsilon=0.1, degree=1)(x))"
" # Call the RBF interpolator\n",
" val = RBFInterpolator(snapshots[:,0:2:1], snapshots[:,2], kernel='multiquadric', epsilon=3)(x)\n",
"\n",
" return val"
]
},
{
Expand All @@ -638,8 +641,8 @@
" x = x.reshape(-1,2) # Reshape array in case its shape is (2,:)\n",
"\n",
" dim = 2\n",
" level = 3\n",
" grid = pysgpp.Grid.createLinearClenshawCurtisBoundaryGrid(dim)\n",
" level = 4\n",
" grid = pysgpp.Grid.createPolyClenshawCurtisBoundaryGrid(dim,degree=2)\n",
" gridStorage = grid.getStorage()\n",
" grid.getGenerator().regular(level)\n",
"\n",
Expand All @@ -655,7 +658,7 @@
" y = np.empty((int(x.size/2), 1))\n",
"\n",
" for n in range(int(x.size/2)):\n",
" p[0] = (x[n,0] - 30.0) / 40.0\n",
" p[0] = (x[n,0] - 0.3) / 0.4\n",
" p[1] = (x[n,1] - 0.3) / 0.5\n",
" opEval = pysgpp.createOperationEvalNaive(grid)\n",
" y[n] = opEval.eval(alpha, p)\n",
Expand Down Expand Up @@ -688,8 +691,9 @@
"xobs = snapshots[:,0:2:1]\n",
"yobs = snapshots[:,2]\n",
"\n",
"xgrid = np.mgrid[30:70:50j, 0.3:0.8:50j]\n",
"xflat = xgrid.reshape(2, -1).T\n",
"xgrid = np.mgrid[0.3:0.7:50j, 0.3:0.8:50j]\n",
"xflat = xgrid.copy()\n",
"xflat = xflat.reshape(2, -1).T\n",
"\n",
"if (metaModel == \"SmolyakInterpolation\"):\n",
" yflat = ResponseSurfaceSmolyak(xflat)\n",
Expand All @@ -699,9 +703,14 @@
"ygrid = yflat.reshape(50, 50)\n",
"\n",
"fig, ax = plt.subplots()\n",
"ax.pcolormesh(*xgrid, ygrid, vmin=200, vmax=1200, shading='gouraud')\n",
"p = ax.scatter(*xobs.T, c=yobs, s=10, ec='k', vmin=200, vmax=1200)\n",
"ax.pcolormesh(*xgrid, ygrid, vmin=350, vmax=450, shading='gouraud')\n",
"p = ax.scatter(*xobs.T, c=yobs, s=10, ec='k', vmin=350, vmax=450)\n",
"\n",
"fig.colorbar(p)\n",
"plt.xlabel('Anastomosis angle')\n",
"plt.ylabel('Residual flow')\n",
"#plt.title('Energy functional E(u)')\n",
"plt.title('Vorticity functional V(u)')\n",
"plt.show()"
]
},
Expand Down Expand Up @@ -781,14 +790,14 @@
},
"outputs": [],
"source": [
"bounds = Bounds([30, 0.3], [70, 0.8])\n",
"bounds = Bounds([0.3, 0.3], [0.7, 0.8])\n",
"\n",
"x0 = np.array([45, 0.4])\n",
"x0 = np.array([0.65, 0.45])\n",
"\n",
"if (metaModel == \"SmolyakInterpolation\"):\n",
" res = minimize(ResponseSurfaceSmolyak, x0, method='trust-constr', jac=\"2-point\", hess=SR1(), constraints=[], options={'verbose': 1}, bounds=bounds)\n",
" res = minimize(ResponseSurfaceSmolyak, x0, method='trust-constr', constraints=[], options={'verbose': 1}, bounds=bounds, tol=1e-6)\n",
"else:\n",
" res = minimize(ResponseSurfaceRBF, x0, method='trust-constr', jac=\"2-point\", hess=SR1(), constraints=[], options={'verbose': 1}, bounds=bounds)\n",
" res = minimize(ResponseSurfaceRBF, x0, method='trust-constr', constraints=[], options={'verbose': 1}, bounds=bounds, tol=1e-6)\n",
"\n",
"print(res.x)"
]
Expand All @@ -799,7 +808,7 @@
"id": "WQ-ifOwVSHJI"
},
"source": [
"Despite needing a few hundred quasi-Newton iterations, an optimal design is found within less than a second. The optimal angle is $\\theta^*=42.5^{\\circ}$ and the optimal residual flow is $\\omega^*=76.4\\%$.\n",
"Despite needing a few hundred quasi-Newton iterations, an optimal design is found within a few seconds. The optimal angle is $\\theta^*=30.0^{\\circ}$ and the optimal residual flow is $\\omega^*=74.3\\%$ for the functional $E$, and $\\theta^*=30.2^{\\circ}$ and the optimal residual flow is $\\omega^*=74.0\\%$ for the functional $V$.\n",
"\n",
"We can also solve just for the optimal angle $\\theta^*$ is the residual flow $\\omega$ is known by modifying the bound constraints:"
]
Expand All @@ -816,14 +825,14 @@
},
"outputs": [],
"source": [
"bounds = Bounds([30, 0.4], [70, 0.4])\n",
"bounds = Bounds([0.3, 0.4], [0.7, 0.4])\n",
"\n",
"x0 = np.array([45, 0.4])\n",
"x0 = np.array([0.6, 0.4])\n",
"\n",
"if (metaModel == \"SmolyakInterpolation\"):\n",
" res = minimize(ResponseSurfaceSmolyak, x0, method='trust-constr', jac=\"2-point\", hess=SR1(), constraints=[], options={'verbose': 1}, bounds=bounds)\n",
" res = minimize(ResponseSurfaceSmolyak, x0, method='trust-constr', constraints=[], options={'verbose': 1}, bounds=bounds)\n",
"else:\n",
" res = minimize(ResponseSurfaceRBF, x0, method='trust-constr', jac=\"2-point\", hess=SR1(), constraints=[], options={'verbose': 1}, bounds=bounds)\n",
" res = minimize(ResponseSurfaceRBF, x0, method='trust-constr', constraints=[], options={'verbose': 1}, bounds=bounds)\n",
"\n",
"print(res.x)"
]
Expand All @@ -834,7 +843,7 @@
"id": "sniXVYEHS4Ey"
},
"source": [
"Thus for a fixed flow residual $\\omega = 40\\%$, the optimal angle is now $\\theta^*=60.2^{\\circ}$.\n",
"Thus for a fixed flow residual $\\omega = 40\\%$, the optimal angle is now $\\theta^*=45.6^{\\circ}$ for the functional $E$, and $\\theta^*=45.9^{\\circ}$ for the functional $V$.\n",
"\n",
"---"
]
Expand Down

0 comments on commit a4874e8

Please sign in to comment.