{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Eigenvalue problem for phase shift\n",
    "\n",
    "[WIP: Master's thesis A. Schöfl]\n",
    "\n",
    "Let us solve the eigenvalue problem of finding $(k,u)$ satisfying \n",
    "\n",
    "$$\\int_\\Omega \\frac{1}{\\mu} \\nabla \\left( e^{i k x} u(x,y)\\right) \\cdot \\nabla \\left( e^{-i k x} v(x,y)\\right) = \\omega^2 \\int_\\Omega \\epsilon \\, u v $$\n",
    "\n",
    "leads to quadratic eigenvalue problem for $ik$:\n",
    "\n",
    "$$ \n",
    "\\underbrace{\\int_\\Omega  \\frac{1}{\\mu}  \\nabla u \\nabla v -\\omega^2 \\int_\\Omega \\epsilon \\, u v }_A\n",
    "+ i k  \\underbrace{\\int_\\Omega \\frac{1}{\\mu}\\left( u \\partial_x v- v\\partial_x u\\right) }_B \n",
    "+ (ik)^2 \\underbrace{\\int_\\Omega \\frac{-1}{\\mu} u v }_C \n",
    "=0$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve import *\n",
    "from ngsolve.webgui import Draw\n",
    "import math\n",
    "import scipy.linalg\n",
    "import numpy as np\n",
    "from ngsolve.eigenvalues import SOAR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nr_eigs = 200\n",
    "omega=Parameter(5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Periodic unit-cell"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from netgen.occ import *\n",
    "from netgen.webgui import Draw as DrawGeo\n",
    "\n",
    "a = 1\n",
    "r = 0.38\n",
    "\n",
    "rect = WorkPlane().RectangleC(a,a).Face()\n",
    "circ = WorkPlane().Circle(0,0, r).Face()\n",
    "# r2 = WorkPlane().Rotate(30).RectangleC(a, a/10).Face()\n",
    "# circ += r2\n",
    "\n",
    "outer = rect-circ\n",
    "inner = rect*circ\n",
    "\n",
    "outer.faces.name = \"outer\"\n",
    "outer.faces.col=(1,1,0)\n",
    "    \n",
    "inner.faces.col=(1,0,0)\n",
    "inner.faces.name=\"inner\"\n",
    "shape = Glue([outer, inner])\n",
    "\n",
    "shape.edges.Max(X).name = \"right\"\n",
    "shape.edges.Max(-X).name = \"left\"\n",
    "shape.edges.Max(Y).name = \"top\"\n",
    "shape.edges.Max(-Y).name = \"bot\"\n",
    "\n",
    "shape.edges.Max(Y).Identify(shape.edges.Min(Y), \"bt\")\n",
    "shape.edges.Max(X).Identify(shape.edges.Min(X), \"lr\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DrawGeo(shape);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh = Mesh(OCCGeometry(shape, dim=2).GenerateMesh(maxh=0.1))\n",
    "mesh.Curve(5)\n",
    "Draw (mesh);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting up finite element system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fes = Compress(Periodic(H1(mesh, complex=True, order=5)))\n",
    "# fes = Periodic(H1(mesh, complex=True, order=4))\n",
    "\n",
    "print (fes.ndof)\n",
    "u,v = fes.TnT()  \n",
    "\n",
    "cf_mu = 1\n",
    "cf_eps = mesh.MaterialCF({\"inner\":9}, default=1)\n",
    "\n",
    "a = BilinearForm( 1/cf_mu*grad(u)*grad(v)*dx-cf_eps*omega**2*u*v*dx )\n",
    "b = BilinearForm( 1/cf_mu*(u*grad(v)[0]-grad(u)[0]*v)*dx )\n",
    "c = BilinearForm( -1/cf_mu*u*v*dx )\n",
    "a1 = BilinearForm( 1/cf_mu*grad(u)*grad(v)*dx )\n",
    "a2 = BilinearForm( cf_eps*u*v*dx )\n",
    "# a = a1 - omega**2 * a2\n",
    "\n",
    "\n",
    "a.Assemble()\n",
    "b.Assemble()\n",
    "c.Assemble()\n",
    "a1.Assemble()\n",
    "a2.Assemble();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "inv = a.mat.Inverse(freedofs=fes.FreeDofs(), inverse=\"sparsecholesky\")\n",
    "M1 = -inv@b.mat\n",
    "M2 = -inv@c.mat\n",
    "Q = SOAR(M1,M2, nr_eigs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(Q)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Small QEP \n",
    "$$\n",
    "(A_m + \\lambda B_m + \\lambda^2 C_m) x = 0\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def SolveProjectedSmall (Am, Bm, Cm):\n",
    "    half = Am.h\n",
    "    Cmi = Cm.I\n",
    "    K = Matrix(2*half, 2*half, complex = fes.is_complex)\n",
    "\n",
    "    K[:half, :half] = -Cmi*Bm\n",
    "    K[:half, half:] = -Cmi*Am\n",
    "    K[half:, :half] = Matrix(half, complex=True).Identity()\n",
    "    K[half:, half:] = 0\n",
    "\n",
    "    Kr = Matrix(2*half)\n",
    "    Kr.A = K.A.real\n",
    "    lam, eig = scipy.linalg.eig(Kr)\n",
    "    vecs = Matrix(eig)[0:len(Q),:] \n",
    "    lam = Vector(lam)\n",
    "\n",
    "    return lam, vecs    \n",
    "\n",
    "def SolveProjected(mata, matb, matc, Q):\n",
    "    Am = InnerProduct(Q, mata*Q, conjugate = True)\n",
    "    Bm = InnerProduct(Q, matb*Q, conjugate = True)\n",
    "    Cm = InnerProduct(Q, matc*Q, conjugate = True)\n",
    "    return SolveProjectedSmall (Am, Bm, Cm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lams, vecs = SolveProjected(a.mat, b.mat, c.mat, Q)\n",
    "Z = (Q * vecs).Evaluate()\n",
    "for vec in Z: vec /= Norm(vec)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gfu = GridFunction(fes)\n",
    "ind = np.absolute(np.asarray(lams)-6j).argmin()\n",
    "gfu.vec.data = Z[ind]\n",
    "Draw (gfu, animate_complex=True, deformation=True, scale=2);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "fig, ax = plt.subplots()\n",
    "ax.set_xlim(-20,20)\n",
    "ax.set_ylim(-20,20)\n",
    "plt.plot([l.real for l in lams],[l.imag for l in lams],'x')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Experiments with reduced basis:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# help (Q)\n",
    "Qred = MultiVector(Q[0], 0)\n",
    "\n",
    "for fs in np.linspace(0,0.8,5):\n",
    "    omega.Set(2*math.pi*fs)\n",
    "    print (\"fs =\", fs)\n",
    "    a.Assemble()\n",
    "    \n",
    "    inv = a.mat.Inverse(freedofs=fes.FreeDofs(), inverse=\"sparsecholesky\")\n",
    "    M1 = -inv@b.mat\n",
    "    M2 = -inv@c.mat\n",
    "    Q = SOAR(M1,M2, nr_eigs)    \n",
    "    lams, vecs = SolveProjected(a.mat, b.mat, c.mat, Q)\n",
    "    Z = (Q * vecs).Evaluate()\n",
    "    for vec in Z: vec /= Norm(vec)    \n",
    "\n",
    "    for lam, vec in zip(lams, Z):\n",
    "      if abs(lam) < 6:\n",
    "        if lam.imag > 1e-8:\n",
    "            hvr = vec.CreateVector()\n",
    "            hvi = vec.CreateVector()\n",
    "            for i in range(len(vec)):\n",
    "                hvr[i] = vec[i].real\n",
    "                hvi[i] = vec[i].imag\n",
    "            Qred.Append (hvr)\n",
    "            Qred.Append (hvi)\n",
    "        if abs(lam.imag) < 1e-8:\n",
    "            hvr = vec.CreateVector()\n",
    "            for i in range(len(vec)):\n",
    "                hvr[i] = vec[i].real\n",
    "            Qred.Append (hvr)            \n",
    "\n",
    "Qred.Orthogonalize()\n",
    "lamsred, vecsred = SolveProjected(a.mat, b.mat, c.mat, Qred)\n",
    "\n",
    "print (\"dim Qred=\", len(Qred))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fs = []\n",
    "ks =[]\n",
    "ksi =[]\n",
    "\n",
    "A1m = InnerProduct(Qred, a1.mat*Qred, conjugate = True)\n",
    "A2m = InnerProduct(Qred, a2.mat*Qred, conjugate = True)\n",
    "Bm = InnerProduct(Qred, b.mat*Qred, conjugate = True)\n",
    "Cm = InnerProduct(Qred, c.mat*Qred, conjugate = True)\n",
    "\n",
    "for fi in np.linspace(0, 0.7, 1000): \n",
    "    # print (\"fi =\", fi)\n",
    "    omega.Set(2*math.pi*fi)\n",
    "    Am = A1m - omega.Get()**2 * A2m\n",
    "    lamsred, vecsred = SolveProjectedSmall(Am, Bm, Cm)\n",
    "    \n",
    "    for lamred in lamsred:\n",
    "        if abs(lamred.real) < 2 and lamred.imag >= 0 and lamred.imag < 6.29:\n",
    "            fs.append(fi)\n",
    "            ks.append(lamred.imag)\n",
    "            ksi.append(lamred.real)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Agrees very well with Scheiber et al, path $\\Gamma-X$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.set_xlim(0,6.28)\n",
    "ax.set_ylim(0,0.7)\n",
    "plt.plot(ks, fs, \"*\", ms=2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.set_xlim(0,2)\n",
    "ax.set_ylim(0,0.7)\n",
    "plt.plot(ksi, fs, \"*\", ms=2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "fig = plt.figure(figsize=(10,10))\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "ax.set_zlim(0,0.7)\n",
    "plt.plot(ks, ksi, fs, \"*\", ms=1)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.13.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}