{ "cells": [ { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# 3.4 Nonlinear minimization problems\n", "We consider minimization problems of the form\n", "\n", "$$\\text{find } u \\in V \\text{ s.t. } E(u) \\leq E(v) \\quad \\forall~ v \\in V.$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "code_folding": [ 0 ], "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "# imports\n", "from ngsolve import *\n", "from ngsolve.webgui import Draw" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "-" } }, "source": [ "Similar to the previous unit we want solve the (minimization) problem using Newton's method. However this time we don't start with an equation but a minimization problem. We will let `NGSolve` derive the corresponding expressions for minimization conditions." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "To solve the problem we use the `Variation` integrator of `NGSolve` and formulate the problem through the symbolic description of an energy functional. \n", "\n", "Let $E(u)$ be the energy that is to be minimized for the unknown state $u$. \n", "\n", "Then a necessary optimality condition is that the derivative at the minimizer $u$ in all directions $v$ vanishes, i.e. \n", "$$\n", " \\delta E(u) (v) = 0 \\quad \\forall v \\in V\n", "$$\n", "\n", "At this point we are back in business of the previous unit as we now have an equation that we need to solve." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "-" } }, "source": [ "Let's continue with a concrete example." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "## Scalar minimization problems\n", "As a first example we take $V = H^1_0$ and \n", "\n", "$$\n", "E(u) = \\int_{\\Omega} \\frac12 \\vert \\nabla u \\vert^2 + \\frac{1}{12} u^4 - fu ~dx.\n", "$$\n", "\n", "The minimization is equivalent (due to convexity) to solving the nonlinear PDE (cf. [unit 3.3](../unit-3.3-nonlinear/nonlinear.ipynb) with $f=10$)\n", "\n", "$$\n", " \\delta E(u)(v) = 0 ~ \\forall v \\in V \\quad \\Leftrightarrow \\quad - \\Delta u + \\frac13 u^3 = f \\text{ in } V'\n", "$$\n", "\n", "with $\\Omega = (0,1)^2$." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "code_folding": [ 0 ], "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "# define geometry and generate mesh\n", "from ngsolve import *\n", "from ngsolve.webgui import *\n", "from netgen.occ import *\n", "shape = Rectangle(1,1).Face()\n", "shape.edges.Min(X).name=\"left\"\n", "shape.edges.Max(X).name=\"right\"\n", "shape.edges.Min(Y).name=\"bottom\"\n", "shape.edges.Max(Y).name=\"top\"\n", "geom = OCCGeometry(shape, dim=2)\n", "mesh = Mesh(geom.GenerateMesh(maxh=0.3))" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "-" } }, "source": [ "We solve the PDE with a Newton iteration." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "V = H1(mesh, order=4, dirichlet=\".*\")\n", "u = V.TrialFunction()" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "We define the semi-linear form expression through the energy functional using the `Variation`-keyword:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "a = BilinearForm (V, symmetric=True)\n", "a += Variation ( (0.5*grad(u)*grad(u) + 1/12*u**4-10*u) * dx)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Now `NGSolve` applies the derivative on the functional so that the previous statement corresponds to:\n", "```\n", "a += (grad(u) * grad(v) + 1/3*u**3*v - 10 * v)*dx\n", "``` \n", "(which has the same form as the problems in [the nonlinear example](../unit-3.7-nonlinear/nonlinear.ipynb))" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "We recall the Newton iteration (cf. [unit-3.3](../unit-3.3-nonlinear/nonlinear.ipynb) ) and now apply the same loop:\n", "\n", "- Given an initial guess $u^0$\n", "- loop over $i=0,..$ until convergence:\n", " - Compute linearization: $A (u^i) + \\delta A(u^i) \\Delta u^{i} = 0 $:\n", " - $f^i = A (u^i)$ \n", " - $B^i = \\delta A(u^i)$ \n", " - Solve $B^i \\Delta u^i = -f^i$\n", " - Update $u^{i+1} = u^i + \\Delta u^{i}$\n", " - Evaluate stopping criteria\n", "- Evaluate $E(u^{i+1})$\n", "\n", "As a stopping criteria we take $\\langle A u^i,\\Delta u^i \\rangle = \\langle A u^i, A u^i \\rangle_{(B^i)^{-1}}< \\varepsilon$." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Note that $A(u^i)$ (`a.Apply(...)`) and $\\delta A(u^i)$ (`a.AssembleLinearization(...)`) are now derived from $A$ which is defined as $A = \\delta E$ through the energy functional $E(\\cdot)$. \n", "\n", "We obtain a similar Newton solver with the two additional advantages:\n", "\n", " * We don't have to form $\\delta E$ manually, but let `NGSolve` do the job and\n", " * we can use the energy functional to interprete the success of iteration steps " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "code_folding": [], "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "def SolveNonlinearMinProblem(a,gfu,tol=1e-13,maxits=10, callback=lambda gfu: None):\n", " res = gfu.vec.CreateVector()\n", " du = gfu.vec.CreateVector()\n", " callback(gfu)\n", " for it in range(maxits):\n", " print (\"Newton iteration {:3}\".format(it),\n", " \", energy = {:16}\".format(a.Energy(gfu.vec)),end=\"\")\n", " \n", " #solve linearized problem:\n", " a.Apply (gfu.vec, res)\n", " a.AssembleLinearization (gfu.vec)\n", " du.data = a.mat.Inverse(V.FreeDofs()) * res\n", " \n", " #update iteration\n", " gfu.vec.data -= du\n", " callback(gfu)\n", "\n", " #stopping criteria\n", " stopcritval = sqrt(abs(InnerProduct(du,res)))\n", " print (\"_{-1}^0.5 = \", stopcritval)\n", " if stopcritval < tol:\n", " break\n", " Redraw(blocking=True)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "So, let's try it out:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "gfu = GridFunction(V)\n", "gfu.Set((x*(1-x))**4*(y*(1-y))**4) # initial guess\n", "gfu_it = GridFunction(gfu.space,multidim=0)\n", "cb = lambda gfu : gfu_it.AddMultiDimComponent(gfu.vec) # store current state\n", "SolveNonlinearMinProblem(a,gfu, callback = cb)\n", "print (\"energy = \", a.Energy(gfu.vec)) " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "Draw(gfu,mesh,\"u\", deformation = True)\n", "#Draw(gfu_it,mesh,\"u\", deformation = True)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Again, a Newton for minimization is also shipped with NGSolve (actually it is the same with the additional knowledge about the Energy and the possibility to do a line search for a given search direction):" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "from ngsolve.solvers import *\n", "gfu.Set((x*(1-x))**4*(y*(1-y))**4) # initial guess\n", "NewtonMinimization(a,gfu)\n", "#Draw(gfu,mesh,\"u\", deformation = True)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Nonlinear elasticity\n", "\n", "We consider a beam which is fixed on one side and is subject to gravity only. We assume a Neo-Hookean hyperelastic material. The model is a nonlinear minimization problem with \n", "\n", "$$\n", " E(v) := \\int_{\\Omega} \\frac{\\mu}{2} ( \\operatorname{tr}(F^T F-I)+\\frac{2 \\mu}{\\lambda} \\operatorname{det}(F^T F)^{\\frac{\\lambda}{2\\mu}} - 1) - \\gamma ~ (f,v) ~~ dx\n", "$$\n", "\n", "where $\\mu$ and $\\lambda$ are the Lamé parameters and $F = I + D v$ where $v: \\Omega \\to \\mathbb{R}^2$ is the sought for displacement." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "We fix the domain to $\\Omega = (0,1) \\times (0,0.1)$" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "code_folding": [ 0 ], "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "# define geometry and generate mesh\n", "shape = Rectangle(1,0.1).Face()\n", "shape.edges.Min(X).name=\"left\"\n", "shape.edges.Min(X).maxh=0.01\n", "shape.edges.Max(X).name=\"right\"\n", "shape.edges.Min(Y).name=\"bot\"\n", "shape.edges.Max(Y).name=\"top\"\n", "geom = OCCGeometry(shape, dim=2)\n", "mesh = Mesh(geom.GenerateMesh(maxh=0.05))" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "# E module and poisson number:\n", "E, nu = 210, 0.2\n", "# Lamé constants:\n", "mu = E / 2 / (1+nu)\n", "lam = E * nu / ((1+nu)*(1-2*nu))\n", "\n", "V = VectorH1(mesh, order=2, dirichlet=\"left\")\n", "u = V.TrialFunction()\n", "\n", "#gravity:\n", "force = CoefficientFunction( (0,-1) )" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Now, we recall the energy\n", "$$\n", " E(v) := \\int_{\\Omega} \\frac{\\mu}{2} ( \\operatorname{tr}(F^T F-I)+\\frac{2 \\mu}{\\lambda} \\operatorname{det}(F^T F)^{\\frac{\\lambda}{2\\mu}} - 1) - \\gamma ~ (f,v) ~~ dx\n", "$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "def Pow(a, b):\n", " return exp (log(a)*b)\n", " \n", "def NeoHook (C):\n", " return 0.5 * mu * (Trace(C-I) + 2*mu/lam * Pow(Det(C), -lam/2/mu) - 1)\n", "\n", "I = Id(mesh.dim)\n", "F = I + Grad(u)\n", "C = F.trans * F\n", "\n", "factor = Parameter(1.0)\n", "\n", "a = BilinearForm(V, symmetric=True)\n", "a += Variation( NeoHook (C).Compile() * dx \n", " -factor * (InnerProduct(force,u) ).Compile() * dx)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "We want to solve the minimization problem for $\\gamma = 5$. Due to the high nonlinearity in the problem, the Newton iteration will not convergence for any initial guess. We approach the case $\\gamma = 5$ by solving problems with $\\gamma = i/10$ for $i=1,..,50$ and taking the solution of the previous problem as an initial guess." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "gfu = GridFunction(V)\n", "gfu.vec[:] = 0\n", "gfu_l = GridFunction(V,multidim=0)\n", "gfu_l.AddMultiDimComponent(gfu.vec)\n", "for loadstep in range(50):\n", " print (\"loadstep\", loadstep)\n", " factor.Set ((loadstep+1)/10)\n", " SolveNonlinearMinProblem(a,gfu)\n", " if (loadstep + 1) % 10 == 0:\n", " gfu_l.AddMultiDimComponent(gfu.vec)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "Draw(gfu_l,mesh, interpolate_multidim=True, animate=True, \n", " deformation=True, min=0, max=1, autoscale=False)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Supplementary 1: Allen-Cahn equation\n", "\n", "The Allen-Cahn equations describe the process of phase separation and is the ($L^2$) gradient-flow equation to the energy\n", "$$\n", " E(v) = \\int_{\\Omega} \\varepsilon \\vert \\nabla v \\vert^2~+~ \\underbrace{v^2(1-v^2)}_{\\Psi(v)} ~ dx\n", "$$\n", "where $\\Psi(v)$ is the so-called double-well potential with the minima $-1,0,1$ (with $0$ being an unstable minima).\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "\n", "The solution to the Allen-Cahn equation solves\n", "\n", "$$\n", "\\partial_t u = \\frac{\\delta E}{\\delta u}\n", "$$\n", "\n", "The quantity $u$ is an indicator for a phase where $-1$ refers to one phase and $1$ to another phase. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "The equation has two driving forces:\n", "\n", "- $u$ is pulled into one of the two stable minima states ($-1$ and $1$) of the nonlinear term $u^2(1-u^2)$ (separation of the phases)\n", "- the diffusion term scaled with $\\varepsilon$ enforces a smooth transition between the two phases. $\\varepsilon$ determines the size of the transition layer\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "-" } }, "source": [ "\n", "We use the `Energy` formulation for energy minimization combined with a simple time stepping with an implicit Euler discretization:\n", "\n", "$$\n", " M u^{n+1} - M u^n = \\Delta t \\underbrace{\\frac{\\delta E}{\\delta u}}_{=:A(u)} (u^{n+1})\n", "$$\n", "\n", "which we can interprete as a nonlinear minimization problem again with the energy\n", "\n", "$$\n", " E^{IE}(v) = \\int_{\\Omega} \\frac{\\varepsilon}{2} \\vert \\nabla v \\vert^2~+~v^2(1-v^2) + \\frac{1}{2\\Delta t} \\vert v - u^n \\vert^2 ~ dx\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "\n", "To solve the nonlinear equation at every time step we again rely on Newton's method.\n", "We first define the periodic geometry, setup the formulation and then apply Newton's method in the next cells:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "# define periodic geometry and generate mesh\n", "shape = Rectangle(1,1).Face()\n", "right=shape.edges.Max(X)\n", "right.name=\"right\"\n", "shape.edges.Min(X).Identify(right,name=\"left\")\n", "top=shape.edges.Max(Y)\n", "top.name=\"top\"\n", "shape.edges.Min(Y).Identify(top,name=\"bottom\")\n", "geom = OCCGeometry(shape, dim=2)\n", "mesh = Mesh(geom.GenerateMesh(maxh=0.1))\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "#use a periodic fe space correspondingly\n", "V = Periodic(H1(mesh, order=3))\n", "u = V.TrialFunction()\n", "\n", "eps = 4e-3\n", "dt = Parameter(1e-1)\n", "gfu = GridFunction(V)\n", "gfuold = GridFunction(V)\n", "a = BilinearForm (V, symmetric=False)\n", "a += Variation( (eps/2*grad(u)*grad(u) + ((1-u*u)*(1-u*u)) \n", " + 0.5/dt*(u-gfuold)*(u-gfuold)) * dx)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "from math import pi\n", "gfu = GridFunction(V)\n", "#gfu.Set(sin(2*pi*x)) # regular initial values\n", "gfu.Set(sin(1e7*(x+y*y))) #<- essentially a random function\n", "gfu_t = GridFunction(V, multidim=0)\n", "gfu_t.AddMultiDimComponent(0.1*gfu.vec)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "t = 0; tend = 5; cnt = 0; sample_rate = int(floor(0.5/dt.Get()))\n", "while t < tend - 0.5 * dt.Get():\n", " gfuold.vec.data = gfu.vec\n", " SolveNonlinearMinProblem(a,gfu)\n", " if (cnt+1) % sample_rate == 0:\n", " gfu_t.AddMultiDimComponent(0.1*gfu.vec)\n", " t += dt.Get(); cnt += 1\n", " print(\"t = \", t)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "Draw(gfu_t, mesh, interpolate_multidim=True, animate=True,\n", " min=-0.1, max=0.1, autoscale=False, deformation=True)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Supplementary 2: Minimal energy extension (postscript in [unit-2.1.3](../unit-2.1.3-bddc/bddc.ipynb) )\n", "\n", "In [unit-2.1.3](../unit-2.1.3-bddc/bddc.ipynb) we discussed the BDDC preconditioner and characterized the coarse grid solution as the condensed problem with the continuity only w.r.t. the coarse grid dofs.\n", "\n", "We can characterize this as a minimization problem involving\n", "\n", "$$\n", " u \\in V^{ho,disc}, \\quad u^{lo,cont} \\in V^{lo,cont}, \\quad \\lambda \\in V^{lo,disc},\n", "$$\n", "\n", "Here $u$ is the solution to the coarse-grid (decoupled) problem where $u^{lo,cont}$ represents the coarse (continuity-)dofs and $\\lambda$ is the corresponding Lagrange multiplier to enforce continuity of $u$ w.r.t. the `lo,cont`-dofs.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "mesh = Mesh(geom.GenerateMesh(maxh=0.1))\n", "fes_ho = Discontinuous(H1(mesh, order=10))\n", "fes_lo = H1(mesh, order=1, dirichlet=\".*\")\n", "fes_lam = Discontinuous(H1(mesh, order=1))\n", "fes = fes_ho*fes_lo*fes_lam\n", "uho, ulo, lam = fes.TrialFunction()" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "The energy that is to be minimized is:\n", "\n", "$$\n", "\\int_{\\Omega} \\frac12 \\Vert \\nabla u \\Vert^2 - u + \\sum_T \\sum_{V \\in V(T)} ((u-u^{lo})\\cdot \\lambda)|_{V} \\longrightarrow \\operatorname{min}!\n", "$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [], "source": [ "a = BilinearForm(fes)\n", "a += Variation(0.5 * grad(uho)*grad(uho)*dx \n", " - 1*uho*dx \n", " + (uho-ulo)*lam*dx(element_vb=BBND))\n", "gfu = GridFunction(fes)\n", "NewtonMinimization(a=a, u=gfu)\n", "Draw(gfu.components[0],mesh,deformation=True)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "The minimization problem is solved by the solution of the PDE:\n", "$$\n", "\\int_{\\Omega} \\nabla u \\cdot \\nabla v = \\int_{\\Omega} 1 \\cdot v \\quad \\forall ~ v \\in V^{ho,disc}\n", "$$\n", "under the constraint\n", "$$\n", " u(v) = u^{lo}(v) \\quad \\text{ for all vertices } v \\in V(T) \\text{ for all } T.\n", "$$" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.10.10" }, "varInspector": { "cols": { "lenName": 16, "lenType": 16, "lenVar": 40 }, "kernels_config": { "python": { "delete_cmd_postfix": "", "delete_cmd_prefix": "del ", "library": "var_list.py", "varRefreshCmd": "print(var_dic_list())" }, "r": { "delete_cmd_postfix": ") ", "delete_cmd_prefix": "rm(", "library": "var_list.r", "varRefreshCmd": "cat(var_dic_list()) " } }, "types_to_exclude": [ "module", "function", "builtin_function_or_method", "instance", "_Feature" ], "window_display": false } }, "nbformat": 4, "nbformat_minor": 2 }