NGSolve24

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "4307bd80-260f-47c9-91a9-ad0056c7fc47",
   "metadata": {},
   "source": [
    "# Solver Layers\n",
    "\n",
    "partially work in progress ... "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "747c140a-d7db-4181-a497-a56e4357bcb3",
   "metadata": {},
   "source": [
    "**Current situation:**\n",
    "\n",
    "- We implement a time-stepper\n",
    "- inside the time-stepper we run some Newton-like method to solve a non-linear equation.\n",
    "- Newton's method computes residual and tangent, in uses some direct solver to invert.\n",
    "\n",
    "**Layer approach:**\n",
    "\n",
    "- We create a `Preconditoner` object for the `BilinearForm`\n",
    "- We create a `LinearSolver` object from the `BilienarForm` and the `Preconditioner` object\n",
    "- We create a `NonlinearSolver` from the `BilinearForm` and the `LienarSolver`\n",
    "- We create a `TimerStepper` from the "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d588c82-f530-4f37-8ffc-885e9d1ae275",
   "metadata": {},
   "source": [
    "## Example:\n",
    "\n",
    "setting up the model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "7c636c32-feca-436d-b07e-47dfa2fb0f68",
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve import *\n",
    "from ngsolve.webgui import Draw\n",
    "\n",
    "from netgen.occ import *\n",
    "shape = Rectangle(1,0.1).Face()\n",
    "shape.edges.Max(X).name=\"right\"\n",
    "shape.edges.Min(X).name=\"left\"\n",
    "shape.edges.Max(Y).name=\"top\"\n",
    "shape.edges.Min(Y).name=\"bot\"\n",
    "shape.vertices.Min(X+Y).maxh=0.01\n",
    "shape.vertices.Min(X-Y).maxh=0.01\n",
    "mesh = Mesh(OCCGeometry(shape, dim=2).GenerateMesh(maxh=0.05))\n",
    "\n",
    "E, nu = 210, 0.2\n",
    "mu  = E / 2 / (1+nu)\n",
    "lam = E * nu / ((1+nu)*(1-2*nu))\n",
    "\n",
    "def C(u):\n",
    "    F = Id(2) + Grad(u)\n",
    "    return F.trans * F\n",
    "\n",
    "def NeoHooke (C):\n",
    "    return 0.5*mu*(Trace(C-Id(2)) + 2*mu/lam*Det(C)**(-lam/2/mu)-1)\n",
    "\n",
    "factor = Parameter(0.5)\n",
    "force = CoefficientFunction( (0,factor) )\n",
    "\n",
    "fes = H1(mesh, order=4, dirichlet=\"left\", dim=mesh.dim)\n",
    "u,v  = fes.TnT()\n",
    "\n",
    "a = BilinearForm(fes, symmetric=True)\n",
    "a += Variation(NeoHooke(C(u)).Compile()*dx)\n",
    "a += Variation((-InnerProduct(force,u)).Compile()*dx)\n",
    "gfu = GridFunction(fes)\n",
    "gfu.vec[:] = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d56bf1a1-bd50-4ab3-a070-f3b781c183df",
   "metadata": {},
   "source": [
    "### Usual approach\n",
    "\n",
    "Call Newton's method, with a couple of arguments:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "76215d03-9efa-4641-b376-0104e219b70b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c575f8dcad264ecc91d3444843fd03bf",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(layout=Layout(height='50vh', width='100%'), value={'gui_settings': {}, 'ngsolve_version': '6.2.24…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "solvers.Newton(a=a, u=gfu, inverse=\"sparsecholesky\", printing=False)\n",
    "scene = Draw (C(gfu)[0,0]-1, mesh, deformation=gfu, min=-0.1, max=0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "87d1804e-9ad4-4181-9b37-ac21e31502e5",
   "metadata": {},
   "source": [
    "### Solver layers approach:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "c7729fee-7b16-4023-8748-b7fb107ae1cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "pre = preconditioners.MultiGrid(a)\n",
    "linsolve = solvers.CGSolver(a.mat, pre, maxiter=20)\n",
    "nlsolve = nonlinearsolvers.NewtonSolver(a=a, u=gfu, solver=linsolve)\n",
    "\n",
    "gfu.vec[:] = 0\n",
    "nlsolve.Solve(printing=False);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "11ee2a9f-ce6c-4c29-9b62-c942a37034a4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "21ef819b24ad407aba5450a53732bd4d",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(layout=Layout(height='50vh', width='100%'), value={'gui_settings': {}, 'ngsolve_version': '6.2.24…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "scene = Draw (C(gfu)[0,0]-1, mesh, deformation=gfu, min=-0.1, max=0.1);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5620ca48-58f2-4f36-8385-6065ebe114a5",
   "metadata": {},
   "source": [
    "`Variation` stores the expression for the energy term. We can also symbolically differentiate the energy: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "d0d45926-1307-4ce4-9883-9ac7ea2f027c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "99217842135d435489a8c7084349ccc2",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(layout=Layout(height='50vh', width='100%'), value={'gui_settings': {}, 'ngsolve_version': '6.2.24…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "a = BilinearForm(NeoHooke(C(u)).Diff(u, v)*dx - force*v*dx)\n",
    "\n",
    "gfu = GridFunction(fes)\n",
    "\n",
    "gfu.vec[:] = 0\n",
    "a.Assemble()  # we need that to construct a CGSolver object outside (TODO: solver_cls, solver_args)\n",
    "\n",
    "pre = preconditioners.MultiGrid(a)\n",
    "linsolve = solvers.CGSolver(a.mat, pre, maxiter=20)\n",
    "nlsolve = nonlinearsolvers.NewtonSolver(a=a, u=gfu, solver=linsolve)\n",
    "\n",
    "gfu.vec[:] = 0\n",
    "nlsolve.Solve(printing=False);\n",
    "\n",
    "Draw (C(gfu)[0,0]-1, mesh, deformation=gfu, min=-0.1, max=0.1);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "37f39030-1749-4e60-b9bd-7a6628284a86",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "3d81acca3b2548cdafc33927b10649b5",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "WebGuiWidget(layout=Layout(height='50vh', width='100%'), value={'gui_settings': {}, 'ngsolve_version': '6.2.24…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<ngsolve.comp.GridFunction at 0x10741ea50>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "force = CF( (0,-1))\n",
    "eq = NeoHooke(C(u)).Diff(u,v) * dx - force * v * dx + u.dt.dt * v * dx\n",
    "\n",
    "ts = timestepping.Newmark(equation=eq, dt=10e-2)\n",
    "\n",
    "gfu = GridFunction(fes)\n",
    "scene = Draw(gfu, deformation=True)\n",
    "\n",
    "def callback(t, gfu):\n",
    "    # print(\"t = \", t, \"displacement = \", gfu(mesh(1,0.05)))\n",
    "    scene.Redraw()\n",
    "ts.Integrate(gfu, start_time=0, end_time=20, callback=callback)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38850e07-5389-4485-8706-cb1d040cc779",
   "metadata": {},
   "source": [
    "## Timesteppers\n",
    "\n",
    "experiments ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a5b90e8-f0e2-4d49-92ce-64efb7212756",
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))\n",
    "fes = H1(mesh, order=3, dirichlet=\".*\")\n",
    "u,v = fes.TnT()\n",
    "\n",
    "equ = u.dt*v*dx + grad(u)*grad(v)*dx + v*dx"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e63f56fe-c6c2-4349-9a42-426f70459818",
   "metadata": {},
   "source": [
    "this should go into the time-stepping module:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10fe844c-b7fa-4084-b151-738d46144e3a",
   "metadata": {},
   "source": [
    "### implicit Euler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19f12b5d-727a-4757-8760-690464fa5760",
   "metadata": {},
   "outputs": [],
   "source": [
    "allproxies = equ.GetProxies(trial=True)\n",
    "print (\"proxies:\")\n",
    "for proxy in allproxies: print(proxy, \" python dtorder = \", proxy.dt_order)\n",
    "\n",
    "proxiesdt0 = [ u, grad(u) ]\n",
    "proxiesdt1 = u.dt\n",
    "\n",
    "massterm = equ.Replace( { u : 0 , grad(u) : CF((0,0)) } )\n",
    "print (\"mass-term: \", massterm)\n",
    "\n",
    "stiffnessterm = equ.Replace( { u.dt : 0 } )\n",
    "print (\"stiffness-term:\", stiffnessterm)\n",
    "\n",
    "tau = 1e-2\n",
    "gfold = GridFunction(fes)\n",
    "gfu = GridFunction(fes)\n",
    "\n",
    "mstarterm = equ.Replace( { u.dt : 1/tau*(u-gfold) } )\n",
    "print (\"mstar = \", mstarterm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96df0323-64b2-4438-9cd7-6b47d27d72cd",
   "metadata": {},
   "outputs": [],
   "source": [
    "bfmstar = BilinearForm(mstarterm)\n",
    "bfmstar.Assemble()\n",
    "pre = preconditioners.MultiGrid(bfmstar)\n",
    "linsolve = krylovspace.CGSolver(bfmstar.mat, pre, maxiter=20)\n",
    "nlsolve = nonlinearsolvers.NewtonSolver(a=bfmstar, u=gfu, solver=linsolve)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6b7f7a20-717d-4a7a-9b43-dd1cb760533b",
   "metadata": {},
   "outputs": [],
   "source": [
    "scene = Draw(gfu)\n",
    "\n",
    "for i in range(20):\n",
    "    nlsolve.Solve()\n",
    "    gfold.vec[:] = gfu.vec\n",
    "    scene.Redraw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dcc56445-3426-4894-848b-e6479d64b942",
   "metadata": {},
   "outputs": [],
   "source": [
    "time = Parameter(0)\n",
    "ts = timestepping.ImplicitEuler(equation=equ, time=time, dt=1e-2)\n",
    "gfu = GridFunction(fes)\n",
    "scene = Draw(gfu, deformation=True, scale=-10)\n",
    "def callback(t, gfu):\n",
    "    # print(\"t = \", t, \"norm =\", Norm(gfu.vec))\n",
    "    scene.Redraw()\n",
    "ts.Integrate(gfu, start_time=0, end_time=0.4, callback=callback)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "80f8071e-6598-4b98-bcb4-7acde581a571",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}