saddlepoint/fewconstraints.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "rising-platinum",
   "metadata": {},
   "source": [
    "A Small Number of Constraints\n",
    "---\n",
    "\n",
    "There are many examples where we have only a small number of constraints. An example is a three-phase transformer, where the voltage is prescribed for each coil. This gives three scalar constraints to the field problem.\n",
    "\n",
    "We consider now saddle-point problems\n",
    "\n",
    "$$\n",
    "\\left( \\begin{array}{cc}\n",
    "A & B^T \\\\\n",
    "B & \n",
    "\\end{array} \\right) \n",
    "\\left( \\begin{array}{c}\n",
    "u \\\\ p \n",
    "\\end{array} \\right) =\n",
    "\\left( \\begin{array}{c}\n",
    "f \\\\ g \n",
    "\\end{array} \\right), \n",
    "$$\n",
    "where $B \\in {\\mathbb R}^{m \\times n}$ with $m$ small.\n",
    "\n",
    "Let us reduce the problem to the Schur complement system\n",
    "\n",
    "$$\n",
    "S p = B A^{-1} f - g.\n",
    "$$\n",
    "\n",
    "Now, the Schur complement\n",
    "\n",
    "$$\n",
    "S = B A^{-1} B^T \n",
    "$$\n",
    "\n",
    "is a small matrix, which can be actually computed by solving $m$ problems with the row vectors of $B$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "typical-pickup",
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve import *\n",
    "from ngsolve.webgui import Draw\n",
    "from netgen.geom2d import CSG2d, Circle, Rectangle\n",
    "geo = CSG2d()\n",
    "rect = Rectangle( pmin=(0,0), pmax=(10,1), bottom=\"free\", top=\"free\", left=\"left\", right=\"right\" )\n",
    "geo.Add(rect)\n",
    "mesh = Mesh(geo.GenerateMesh(maxh=0.3))\n",
    "Draw(mesh);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "suspended-belly",
   "metadata": {},
   "source": [
    "We consider an elasticity problem on a rectangular domain, fixed at the left end, and loaded by a force at the right boundary:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "compact-attribute",
   "metadata": {},
   "outputs": [],
   "source": [
    "fes = VectorH1(mesh, order=3, dirichlet=\"left\")\n",
    "u,v = fes.TnT()\n",
    "a = BilinearForm(fes)\n",
    "a += InnerProduct(Sym(Grad(u)), Sym(Grad(v)))*dx + div(u)*div(v)*dx\n",
    "a.Assemble()\n",
    "\n",
    "f = LinearForm(1e-3*v[1]*ds(\"right\")).Assemble()\n",
    "gfu = GridFunction(fes)\n",
    "gfu.vec.data = a.mat.Inverse(fes.FreeDofs()) * f.vec\n",
    "Draw (gfu, deformation=True);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "quality-hayes",
   "metadata": {},
   "source": [
    "Now, instead of the loading, we want to prescribe an average displacement at the right end:\n",
    "\n",
    "$$\n",
    "\\frac{1}{|\\Gamma_r|} \\int_{\\Gamma_r} u_x = U_x \\qquad \\text{and} \\qquad\n",
    "\\frac{1}{|\\Gamma_r|} \\int_{\\Gamma_r} u_y = U_y\n",
    "$$\n",
    "\n",
    "these constraints can be described by linear-forms\n",
    "\n",
    "$$\n",
    "b_x(u) = U_x \\qquad \\text{and} \\qquad b_y(u) = U_y\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "indoor-suspension",
   "metadata": {},
   "source": [
    "Assemble the two constraint-vectors, and write into a MultiVector:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "improving-complex",
   "metadata": {},
   "outputs": [],
   "source": [
    "bx = LinearForm(v[0]*ds(\"right\")).Assemble()\n",
    "by = LinearForm(v[1]*ds(\"right\")).Assemble()\n",
    "\n",
    "B = MultiVector(bx.vec, 2)\n",
    "B[0] = bx.vec\n",
    "B[1] = by.vec"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "established-albuquerque",
   "metadata": {},
   "source": [
    "Compute the Schur complement:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "nonprofit-morocco",
   "metadata": {},
   "outputs": [],
   "source": [
    "S = InnerProduct(B, a.mat.Inverse(fes.FreeDofs()) * B)\n",
    "print (S)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "matched-estonia",
   "metadata": {},
   "source": [
    "Solve for the Lagrange parameter. It has the physical meaning of a force:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "executed-aside",
   "metadata": {},
   "outputs": [],
   "source": [
    "Uxy = Vector( (1, 3) )\n",
    "force = S.I * Uxy\n",
    "print (force)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "tutorial-friendship",
   "metadata": {},
   "source": [
    "and reconstruct the solution $u$ from the computed Lagrange parameter:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "imposed-plane",
   "metadata": {},
   "outputs": [],
   "source": [
    "fvec = B * force\n",
    "gfu.vec.data = a.mat.Inverse(fes.FreeDofs()) * fvec\n",
    "Draw (gfu, deformation=True);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "flying-width",
   "metadata": {},
   "source": [
    "This method requires to solve one linear system with right hand side $f$, and $m$ linerar systems for the Schur complement, one for each scalar constraint. We can use a preconditioner for $A$, and solve these $m+1$ systems by conjugate gradient iterations."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "similar-sacramento",
   "metadata": {},
   "source": [
    "Projected preconditioner\n",
    "---\n",
    "\n",
    "Next we present a method to solve these problem at the cost of one unconstrained problem, plus $m$ applications of the preconditioner. \n",
    "\n",
    "It consists of two steps:\n",
    "- Reduction to homogeneous constraints\n",
    "- Construct a preconditioner on the kernel of $B$\n",
    "\n",
    "**Reduction to homogeneous constraints**\n",
    "\n",
    "We solve the problem in two steps: First we find an $u_0$ such that $B u_0 = g$. We could do so by solving\n",
    "\n",
    "$$\n",
    "\\left( \\begin{array}{cc}\n",
    "A & B^T \\\\\n",
    "B & \n",
    "\\end{array} \\right) \n",
    "\\left( \\begin{array}{c}\n",
    "u_0 \\\\ p \n",
    "\\end{array} \\right) =\n",
    "\\left( \\begin{array}{c}\n",
    "0 \\\\ g \n",
    "\\end{array} \\right)\n",
    "$$\n",
    "\n",
    "using the Schur complement. However, this still requires the solution of $m$ problems for the Schur complement.\n",
    "\n",
    "A cheaper way is to solve \n",
    "\n",
    "$$\n",
    "\\left( \\begin{array}{cc}\n",
    "\\hat{A} & B^T \\\\\n",
    "B & \n",
    "\\end{array} \\right) \n",
    "\\left( \\begin{array}{c}\n",
    "u_0 \\\\ p \n",
    "\\end{array} \\right) =\n",
    "\\left( \\begin{array}{c}\n",
    "0 \\\\ g \n",
    "\\end{array} \\right), \n",
    "$$\n",
    "\n",
    "where $\\hat A$ is a preconditioner to $A$. This can be solved by the inexact Schur complement, where $A^{-1}$ is replace by $\\hat A^{-1}$. We use the very cheap Jacobi preconditioner:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "available-duplicate",
   "metadata": {},
   "outputs": [],
   "source": [
    "prea = a.mat.CreateSmoother(fes.FreeDofs())\n",
    "Shat = InnerProduct(B, prea * B)\n",
    "print (Shat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "invisible-wales",
   "metadata": {},
   "outputs": [],
   "source": [
    "force = Shat.I * Uxy\n",
    "\n",
    "fvec = B * force\n",
    "gfu0 = GridFunction(fes)\n",
    "gfu0.vec.data = prea * fvec\n",
    "Draw (gfu0, deformation=True);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "small-sheep",
   "metadata": {},
   "source": [
    "Now we have to solve the correction problem $u = u_0 + u_1$, where $u_1$ solves the saddle-point problem with homogeneous constraint:\n",
    "\n",
    "$$\n",
    "\\left( \\begin{array}{cc}\n",
    "A & B^T \\\\\n",
    "B & \n",
    "\\end{array} \\right) \n",
    "\\left( \\begin{array}{c}\n",
    "u_1 \\\\ p \n",
    "\\end{array} \\right) =\n",
    "\\left( \\begin{array}{c}\n",
    "f - A u_0 \\\\ 0\n",
    "\\end{array} \\right)\n",
    "$$\n",
    "\n",
    "Let $V_0 = \\operatorname{ker} \\{ B \\}$ be the null-space of $B$. Instead of solving the problem with constraints, we could solve the variational problem on the linear sub-space $V_0$: find $u_1 \\in V_0$ such that\n",
    "\n",
    "$$\n",
    "v^T A u_1 = v^T (f - A u_0) \\qquad v \\in V_0\n",
    "$$\n",
    "\n",
    "We solve this problem on the sub-space $V_0$ by a preconditioned conjugate gradient method, where the range of the projected preconditioner-action $\\hat A^{-1}$ spans the sub-space $V_0$:\n",
    "\n",
    "$$\n",
    "Proj(\\hat A^{-1}) : d \\mapsto w \\qquad \\text{such that}\n",
    "$$\n",
    "\n",
    "$$\n",
    "w \\in V_0 :  \\quad v^T \\hat A w = v^T d \\qquad \\forall v \\in V_0\n",
    "$$\n",
    "\n",
    "The projected problem for the preconditioned system can be solved by the saddle point problem\n",
    "\n",
    "$$\n",
    "\\left( \\begin{array}{cc}\n",
    "\\hat{A} & B^T \\\\\n",
    "B & \n",
    "\\end{array} \\right) \n",
    "\\left( \\begin{array}{c}\n",
    "w \\\\ p \n",
    "\\end{array} \\right) =\n",
    "\\left( \\begin{array}{c}\n",
    "d \\\\ 0\n",
    "\\end{array} \\right).\n",
    "$$\n",
    "\n",
    "With the inexact Schur complement $\\hat S = B \\hat A^{-1} B^T$ this is\n",
    "\n",
    "$$\n",
    "\\hat S p = B \\hat A^{-1} d\n",
    "$$\n",
    "\n",
    "and then\n",
    "\n",
    "$$\n",
    "w = \\hat A^{-1} (d - B^T p) = \\hat A^{-1} (I - B^T S^{-1} B \\hat A^{-1})\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "impressive-german",
   "metadata": {},
   "outputs": [],
   "source": [
    "Bmat = BaseMatrix(B)\n",
    "projpre = prea @ (IdentityMatrix() - Bmat @ BaseMatrix(Shat.I) @ Bmat.T @ prea)\n",
    "\n",
    "res = -a.mat * gfu0.vec\n",
    "res = res.Evaluate()\n",
    "\n",
    "from ngsolve.krylovspace import CGSolver\n",
    "projinv = CGSolver(a.mat, projpre, printrates=False, maxiter=1000)\n",
    "\n",
    "gfu1 = GridFunction(fes)\n",
    "gfu1.vec.data = projinv * res\n",
    "Draw (gfu1, deformation=True)\n",
    "\n",
    "gfu1.vec.data += gfu0.vec\n",
    "Draw (gfu1, deformation=True);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "desperate-disclaimer",
   "metadata": {},
   "source": [
    "Some remarks:\n",
    "* instead of applying the preconditioner $\\hat A$ twice, one could store the matrix $\\hat A^{-1} B$\n",
    "\n",
    "* the Lagrange parameter $p$ can be computed by multiplying the equation $A u + B^T p = f$ from left with $B \\hat A^{-1}$, and then solving with the inexact Schur complement\n",
    "$$\n",
    "p = (B \\hat A^{-1} B^T)^{-1} (B \\hat A^{-1} (f - A u))\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "threaded-story",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.9.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}