{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# This file is part of nannos\n# License: GPLv3\n%matplotlib notebook"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Convergence\n\nConvergence of the various FMM formulations.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import time\n\nimport matplotlib.pyplot as plt\n\nimport nannos as nn\n\nbk = nn.backend"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We will study the convergence on a benchmark case from\n:cite:p:`Li1997`.\nFirst we define the main function that performs the simulation.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "wavelength = 1\nsq_size = 1.25 * wavelength\neps_diel = 2.25\n\n\ndef checkerboard_cellA(nh, formulation):\n    d = 2 * sq_size\n    Nx = 2**9\n    Ny = 2**9\n    lattice = nn.Lattice(([d, 0], [0, d]), discretization=(Nx, Ny))\n    pw = nn.PlaneWave(wavelength=wavelength, angles=(0, 0, 0))\n    epsgrid = lattice.ones() * eps_diel\n    sq1 = lattice.square((0.25 * d, 0.25 * d), sq_size)\n    sq2 = lattice.square((0.75 * d, 0.75 * d), sq_size)\n    epsgrid[sq1] = 1\n    epsgrid[sq2] = 1\n\n    sup = lattice.Layer(\"Superstrate\", epsilon=eps_diel)\n    sub = lattice.Layer(\"Substrate\", epsilon=1)\n    st = lattice.Layer(\"Structured\", wavelength)\n    st.epsilon = epsgrid\n\n    sim = nn.Simulation([sup, st, sub], pw, nh, formulation=formulation)\n    # this actually corresponds to order (0,-1) for the other unit cell in [Li1997]\n    order = (-1, -1)\n    R, T = sim.diffraction_efficiencies(orders=True)\n    t = sim.get_order(T, order)\n    return t, sim\n\n\ndef checkerboard_cellB(nh, formulation):\n    d = sq_size * 2**0.5\n    Nx = 2**9\n    Ny = 2**9\n    lattice = nn.Lattice(([d, 0], [0, d]), discretization=(Nx, Ny))\n    pw = nn.PlaneWave(wavelength=wavelength, angles=(0, 45, 0))\n    epsgrid = lattice.ones() * eps_diel\n    sq = lattice.square((0.5 * d, 0.5 * d), sq_size, rotate=45)\n    epsgrid[sq] = 1\n\n    sup = lattice.Layer(\"Superstrate\", epsilon=eps_diel)\n    sub = lattice.Layer(\"Substrate\", epsilon=1)\n    st = lattice.Layer(\"Structured\", wavelength)\n    st.epsilon = epsgrid\n\n    # st.plot()\n    sim = nn.Simulation([sup, st, sub], pw, nh, formulation=formulation)\n    order = (0, -1)\n    R, T = sim.diffraction_efficiencies(orders=True)\n    t = sim.get_order(T, order)\n    return t, sim"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Perform the simulation for different formulations and number\nof retained harmonics:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def plot_cell(sim):\n    axin = plt.gca().inset_axes([0.77, 0.0, 0.25, 0.25])\n    lay = sim.get_layer_by_name(\"Structured\")\n    lay.plot(ax=axin)\n    axin.set_axis_off()\n\n\nNH = [100, 200, 300, 400, 600]\nformulations = [\"original\", \"tangent\", \"pol\", \"jones\"]\n\nfor icell, cell_fun in enumerate([checkerboard_cellA, checkerboard_cellB]):\n    celltype = \"A\" if icell == 0 else \"B\"\n\n    print(\"============================\")\n    print(f\"cell type {celltype}\")\n    print(\"============================\")\n\n    nhs = {f: [] for f in formulations}\n    ts = {f: [] for f in formulations}\n    times = {f: [] for f in formulations}\n\n    for nh in NH:\n\n        print(\"number of harmonics = \", nh)\n\n        for formulation in formulations:\n            t0 = -time.time()\n            t, sim = cell_fun(nh, formulation=formulation)\n            t0 += time.time()\n            print(\"formulation = \", formulation)\n            print(f\"number of harmonics: asked={nh}, actual={sim.nh}\")\n            print(f\"elapsed time = {t0}s\")\n            print(\"T(0,-1) = \", t)\n            print(\"-----------------\")\n            nhs[formulation].append(sim.nh)\n            ts[formulation].append(t)\n            times[formulation].append(t0)\n\n    #########################################################################\n    # Plot the results:\n\n    markers = {\"original\": \"^\", \"tangent\": \"o\", \"jones\": \"s\", \"pol\": \"^\"}\n    colors = {\n        \"original\": \"#d4b533\",\n        \"tangent\": \"#d46333\",\n        \"jones\": \"#3395d4\",\n        \"pol\": \"#54aa71\",\n    }\n\n    plt.figure()\n    for formulation in formulations:\n        plt.plot(\n            nhs[formulation],\n            ts[formulation],\n            \"-\",\n            color=colors[formulation],\n            marker=markers[formulation],\n            label=formulation,\n        )\n        plt.pause(0.1)\n    plt.legend(loc=5, ncols=2)\n    plt.xlabel(\"number of Fourier harmonics $n_h$\")\n    plt.ylabel(\"$T_{0,-1}$\")\n    # plt.ylim(0.1255, 0.129)\n    plt.title(f\"cell {celltype}\")\n    plot_cell(sim)\n    plt.tight_layout()\n\n    plt.figure()\n\n    for formulation in formulations:\n        plt.plot(\n            nhs[formulation],\n            times[formulation],\n            \"-\",\n            color=colors[formulation],\n            marker=markers[formulation],\n            label=formulation,\n        )\n        plt.pause(0.1)\n    plt.yscale(\"log\")\n    plt.legend(ncols=2)\n    plt.xlabel(\"number of Fourier harmonics $n_h$\")\n    plt.ylabel(\"CPU time (s)\")\n    plt.title(f\"cell {celltype}\")\n    plot_cell(sim)\n    plt.tight_layout()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import nannos.utils.jupyter\n%nannos_version_table"
      ]
    }
  ],
  "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.12.5"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}