NGSolve24

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "4307bd80-260f-47c9-91a9-ad0056c7fc47",
   "metadata": {},
   "source": [
    "# Solver Layers"
   ]
  },
  {
   "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, and uses some direct solver to invert.\n",
    "\n",
    "**Layer approach:**\n",
    "\n",
    "* We create a `Preconditoner` object for the `BilinearForm` object\n",
    "* We create a `LinearSolver` object from the `BilinearForm.mat` and the `Preconditioner` object\n",
    "* We create a `NonlinearSolver` from the `BilinearForm` and the `LinearSolver` object (or LinearSolver-class+args)\n",
    "* We create a `TimeStepper` from a `BilinearForm` object and `NonLinearSolver` class+args"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d588c82-f530-4f37-8ffc-885e9d1ae275",
   "metadata": {},
   "source": [
    "## Example:\n",
    "\n",
    "setting up the model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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",
    "shape.vertices.Min(X+Y).hpref=1\n",
    "shape.vertices.Min(X-Y).hpref=1\n",
    "\n",
    "mesh = Mesh(OCCGeometry(shape, dim=2).GenerateMesh(maxh=0.05))\n",
    "mesh.RefineHP(2)\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": null,
   "id": "76215d03-9efa-4641-b376-0104e219b70b",
   "metadata": {},
   "outputs": [],
   "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": null,
   "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": null,
   "id": "11ee2a9f-ce6c-4c29-9b62-c942a37034a4",
   "metadata": {},
   "outputs": [],
   "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` keeps the expression for the energy. This is not yet supported for `timestepping`. We can use symbolic differentiation to convert it to variational forms:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0d45926-1307-4ce4-9883-9ac7ea2f027c",
   "metadata": {},
   "outputs": [],
   "source": [
    "a = BilinearForm(NeoHooke(C(u)).Diff(u, v)*dx - force*v*dx)\n",
    "gfu = GridFunction(fes)\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",
    "# nlsolve = nonlinearsolvers.NewtonSolver(a=a, u=gfu, lin_solver_cls=solvers.CGSolver, lin_solver_args={ \"pre\":pre })\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": "markdown",
   "id": "c7895515-603b-49f6-911a-1f8c619bc669",
   "metadata": {},
   "source": [
    "## Timesteppers\n",
    "\n",
    "* User provides the time-dependent equation as a variational form\n",
    "* Brand new: time-derivative `u.dt`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37f39030-1749-4e60-b9bd-7a6628284a86",
   "metadata": {},
   "outputs": [],
   "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=1e-1)\n",
    "\n",
    "gfu = GridFunction(fes)\n",
    "scene = Draw(gfu, deformation=True, center=(0,-0.3), radius=0.6)\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=10, callback=callback);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38850e07-5389-4485-8706-cb1d040cc779",
   "metadata": {},
   "source": [
    "### New features to work on exptression trees"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15e0458a-0f99-434c-8116-5980dd9f3aaa",
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))\n",
    "fes = H1(mesh, order=3, dirichlet=\".*\")\n",
    "u,v = fes.TnT()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "310e590c-0d83-449e-9a98-e6c32b32f422",
   "metadata": {},
   "outputs": [],
   "source": [
    "equ = u.dt*v*dx + 1e-3*u*v*dx + grad(u)*grad(v)*dx - v*dx\n",
    "print (equ)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5769646b-65e6-4ec1-a792-99b138abe0d6",
   "metadata": {},
   "source": [
    "we can extract all proxies (trial- or test-functions) from the variational form:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a1ce7430-e743-40a2-b733-60e554307c0d",
   "metadata": {},
   "outputs": [],
   "source": [
    "allproxies = equ.GetProxies(trial=True)\n",
    "print (\"proxies:\")\n",
    "for proxy in allproxies: print(proxy, \" python dtorder = \", proxy.dt_order)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f060090c-2de2-4639-aa15-0b60a8bf0581",
   "metadata": {},
   "outputs": [],
   "source": [
    "proxiesdt0 = [ prox for prox in allproxies if prox.dt_order==0 ]\n",
    "proxiesdt1 = [ prox for prox in allproxies if prox.dt_order==1 ]\n",
    "\n",
    "print (\"order 0 proxies:\")\n",
    "print (*proxiesdt0)\n",
    "\n",
    "print (\"order 1 proxies:\")\n",
    "print (*proxiesdt1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "67772d80-96f7-45b7-b75c-558864751c29",
   "metadata": {},
   "source": [
    "The time-dependent equation is\n",
    "\n",
    "$$\n",
    "\\int \\partial_t u v + uv + \\nabla u \\nabla v - 1 v  \\; dx = 0\n",
    "$$\n",
    "\n",
    "An implicit Euler time-stepping method is to solve \n",
    "\n",
    "$$\n",
    "\\int \\frac{u - u_n}{\\tau} v + uv + \\nabla u \\nabla v - 1 v  \\; dx = 0\n",
    "$$\n",
    "\n",
    "where the unknown $u$ is the value of the new time-step $u_{n+1}$.\n",
    "\n",
    "So we have to replace the \n",
    "\n",
    "$$\n",
    "\\partial_t u  \\rightarrow \\frac{u - u_n}{\\tau}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "864f4bee-1dde-412a-8a8b-5fd16afcbcc9",
   "metadata": {},
   "outputs": [],
   "source": [
    "tau = Parameter(0.1)\n",
    "gfuold = GridFunction(fes)\n",
    "\n",
    "repl = {}\n",
    "for prox in allproxies:\n",
    "    if prox.dt_order==1:\n",
    "        repl[prox] = 1/tau*(prox.anti_dt - prox.anti_dt.ReplaceFunction(gfuold))\n",
    "\n",
    "for key,val in repl.items():\n",
    "    print (\"replace:\", key)\n",
    "    print (\"by\\n\", val)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db77e618-910d-485e-ad9f-5ec3e5504998",
   "metadata": {},
   "outputs": [],
   "source": [
    "ImplEuler_equ = equ.Replace(repl)\n",
    "print (ImplEuler_equ)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f445a879-6ec9-491e-a6bc-b18e0f4df45c",
   "metadata": {},
   "outputs": [],
   "source": [
    "bfIE = BilinearForm(ImplEuler_equ)\n",
    "bfIE.Assemble()\n",
    "pre = preconditioners.MultiGrid(bfIE)\n",
    "linsolve = krylovspace.CGSolver(bfIE.mat, pre, maxiter=20)\n",
    "\n",
    "gfu = GridFunction(fes)\n",
    "nlsolve = nonlinearsolvers.NewtonSolver(a=bfIE, u=gfu, solver=linsolve)\n",
    "\n",
    "scene = Draw(gfu, deformation=True, scale=10)\n",
    "gfuold.vec[:] = 0\n",
    "\n",
    "from time import sleep\n",
    "for i in range(20):\n",
    "    nlsolve.Solve()\n",
    "    gfuold.vec[:] = gfu.vec\n",
    "    scene.Redraw()\n",
    "    sleep(0.2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2c0df60-ebe9-47e9-a367-9ba86cfc1d8d",
   "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
}