{ "cells": [ { "cell_type": "markdown", "id": "696b4bb6", "metadata": {}, "source": [ "# Scattering from a sphere in the axisymmetric formulation\n", "\n", "Copyright (C) 2022 Michele Castriotta, Igor Baratta, Jørgen S. Dokken\n", "\n", "This demo is implemented in three files: one for the mesh\n", "generation with gmsh, one for the calculation of analytical efficiencies,\n", "and one for the variational forms and the solver. It illustrates how to:\n", "\n", "- Setup and solve Maxwell's equations for axisymmetric geometries\n", "- Implement (axisymmetric) perfectly matched layers\n", "\n", "## Equations, problem definition and implementation\n", "\n", "First of all, let's import the modules that will be used:" ] }, { "cell_type": "code", "execution_count": null, "id": "9a30e4d4", "metadata": {}, "outputs": [], "source": [ "import sys\n", "from functools import partial\n", "\n", "import numpy as np\n", "from mesh_sphere_axis import generate_mesh_sphere_axis\n", "from scipy.special import jv, jvp\n", "from dolfinx import fem, mesh, plot\n", "from dolfinx.io import VTXWriter\n", "from dolfinx.io.gmshio import model_to_mesh\n", "from mpi4py import MPI\n", "import ufl\n", "from petsc4py import PETSc\n", "\n", "try:\n", " import gmsh\n", "except ModuleNotFoundError:\n", " print(\"This demo requires gmsh to be installed\")\n", " sys.exit(0)\n", "\n", "try:\n", " import pyvista\n", " have_pyvista = True\n", "except ModuleNotFoundError:\n", " print(\"pyvista and pyvistaqt are required to visualise the solution\")\n", " have_pyvista = False" ] }, { "cell_type": "markdown", "id": "5950420b", "metadata": {}, "source": [ "Since we want to solve time-harmonic Maxwell's equation, we need to solve a\n", "complex-valued PDE, and therefore need to use PETSc compiled with complex numbers." ] }, { "cell_type": "code", "execution_count": null, "id": "8aca707d", "metadata": {}, "outputs": [], "source": [ "if not np.issubdtype(PETSc.ScalarType, np.complexfloating):\n", " print(\"Demo should only be executed with DOLFINx complex mode\")\n", " exit(0)" ] }, { "cell_type": "markdown", "id": "1ffcdf17", "metadata": {}, "source": [ "Now, let's formulate our problem.\n", "Let's consider a metallic sphere immersed in\n", "a background medium (e.g. vacuum or water) and hit by a plane wave.\n", "We want to calculate the scattered electric field scattered.\n", "Even though the problem is three-dimensional,\n", "we can simplify it into many two-dimensional problems\n", "by exploiting its axisymmetric nature. To verify this, let's consider\n", "the weak form of the problem with PML:\n", "\n", "$$\n", "\\begin{align}\n", "&\\int_{\\Omega_m\\cup\\Omega_b}-(\\nabla \\times \\mathbf{E}_s)\n", "\\cdot (\\nabla \\times \\bar{\\mathbf{v}})+\\varepsilon_{r} k_{0}^{2}\n", "\\mathbf{E}_s \\cdot \\bar{\\mathbf{v}}+k_{0}^{2}\\left(\\varepsilon_{r}\n", "-\\varepsilon_b\\right)\\mathbf{E}_b \\cdot \\bar{\\mathbf{v}}~\\mathrm{d} x\\\\\n", "+&\\int_{\\Omega_{pml}}\\left[\\boldsymbol{\\mu}^{-1}_{pml} \\nabla \\times \\mathbf{E}_s\n", "\\right]\\cdot \\nabla \\times \\bar{\\mathbf{v}}-k_{0}^{2}\n", "\\left[\\boldsymbol{\\varepsilon}_{pml} \\mathbf{E}_s \\right]\\cdot\n", "\\bar{\\mathbf{v}}~ d x=0\n", "\\end{align}\n", "$$\n", "\n", "Since we want to exploit the axisymmetry of the problem,\n", "we can use cylindrical coordinates as a more appropriate\n", "reference system:\n", "\n", "$$\n", "\\begin{align}\n", "&\\int_{\\Omega_{cs}}\\int_{0}^{2\\pi}-(\\nabla \\times \\mathbf{E}_s)\n", "\\cdot (\\nabla \\times \\bar{\\mathbf{v}})+\\varepsilon_{r} k_{0}^{2}\n", "\\mathbf{E}_s \\cdot \\bar{\\mathbf{v}}+k_{0}^{2}\\left(\\varepsilon_{r}\n", "-\\varepsilon_b\\right)\\mathbf{E}_b \\cdot\n", "\\bar{\\mathbf{v}}~ \\rho d\\rho dz d \\phi\\\\\n", "+&\\int_{\\Omega_{cs}}\n", "\\int_{0}^{2\\pi}\\left[\\boldsymbol{\\mu}^{-1}_{pml} \\nabla \\times \\mathbf{E}_s\n", "\\right]\\cdot \\nabla \\times \\bar{\\mathbf{v}}-k_{0}^{2}\n", "\\left[\\boldsymbol{\\varepsilon}_{pml} \\mathbf{E}_s \\right]\\cdot\n", "\\bar{\\mathbf{v}}~ \\rho d\\rho dz d \\phi=0\n", "\\end{align}\n", "$$\n", "\n", "Let's now expand $\\mathbf{E}_s$,\n", "$\\mathbf{E}_b$ and $\\bar{\\mathbf{v}}$ in cylindrical harmonics:\n", "\n", "$$\n", "\\begin{align}\n", "& \\mathbf{E}_s(\\rho, z, \\phi) = \\sum_m\\mathbf{E}^{(m)}_s(\\rho, z)e^{-im\\phi} \\\\\n", "& \\mathbf{E}_b(\\rho, z, \\phi) = \\sum_m\\mathbf{E}^{(m)}_b(\\rho, z)e^{-im\\phi} \\\\\n", "& \\bar{\\mathbf{v}}(\\rho, z, \\phi) =\n", "\\sum_m\\bar{\\mathbf{v}}^{(m)}(\\rho, z)e^{+im\\phi}\\\\\n", "\\end{align}\n", "$$\n", "\n", "The curl operator $\\nabla\\times$ in cylindrical coordinates becomes:\n", "\n", "$$\n", "\\begin{aligned}\n", "\\nabla \\times \\mathbf{a}=\n", "\\sum_{m}\\left(\\nabla \\times \\mathbf{a}^{(m)}\\right) e^{-i m \\phi}\n", "\\end{aligned}\n", "$$\n", "\n", "with:\n", "\n", "$$\n", "\\begin{align}\n", "\\left(\\nabla \\times \\mathbf{a}^{(m)}\\right) = &\\left[\\hat{\\rho}\n", "\\left(-\\frac{\\partial a_{\\phi}^{(m)}}{\\partial z}\n", "-i \\frac{m}{\\rho} a_{z}^{(m)}\\right)+\\\\ \\hat{\\phi}\n", "\\left(\\frac{\\partial a_{\\rho}^{(m)}}{\\partial z}\n", "-\\frac{\\partial a_{z}^{(m)}}{\\partial \\rho}\\right)+\\right.\\\\\n", "&\\left.+\\hat{z}\\left(\\frac{a_{\\phi}^{(m)}}{\\rho}\n", "+\\frac{\\partial a_{\\phi}^{(m)}}{\\partial \\rho}\n", "+i \\frac{m}{\\rho} a_{\\rho}^{(m)}\\right)\\right]\n", "\\end{align}\n", "$$\n", "\n", "By implementing these formula in our weak form, and by assuming an\n", "axisymmetric permittivity $\\varepsilon(\\rho, z)$, we can write:\n", "\n", "$$\n", "\\begin{align}\n", "\\sum_{n, m}\\int_{\\Omega_{cs}}&-(\\nabla \\times \\mathbf{E}^{(m)}_s)\n", "\\cdot (\\nabla \\times \\bar{\\mathbf{v}}^{(m)})+\\varepsilon_{r} k_{0}^{2}\n", "\\mathbf{E}^{(m)}_s \\cdot \\bar{\\mathbf{v}}^{(m)}+k_{0}^{2}\n", "\\left(\\varepsilon_{r}\n", "-\\varepsilon_b\\right)\\mathbf{E}^{(m)}_b \\cdot \\bar{\\mathbf{v}}^{(m)}\\\\\n", "&+\\left(\\boldsymbol{\\mu}^{-1}_{pml} \\nabla \\times \\mathbf{E}^{(m)}_s\n", "\\right)\\cdot \\nabla \\times \\bar{\\mathbf{v}}^{(m)}-k_{0}^{2}\n", "\\left(\\boldsymbol{\\varepsilon}_{pml} \\mathbf{E}^{(m)}_s \\right)\\cdot\n", "\\bar{\\mathbf{v}}^{(m)}~ \\rho d\\rho dz \\int_{0}^{2 \\pi} e^{-i(m-n) \\phi}\n", "d \\phi=0\n", "\\end{align}\n", "$$\n", "\n", "For the last integral, we have that:\n", "\n", "$$\n", "\\int_{0}^{2 \\pi} e^{-i(m-n) \\phi}d \\phi=\n", "\\begin{cases}\n", "\\begin{align}\n", "2\\pi ~~~&\\textrm{if } m=n\\\\\n", "0 ~~~&\\textrm{if }m\\neq n\\\\\n", "\\end{align}\n", "\\end{cases}\n", "$$\n", "\n", "We can therefore consider only the case $m = n$ and simplify the summation\n", "in the following way:\n", "\n", "$$\n", "\\begin{align}\n", "\\sum_{m}\\int_{\\Omega_{cs}}&-(\\nabla \\times \\mathbf{E}^{(m)}_s)\n", "\\cdot (\\nabla \\times \\bar{\\mathbf{v}}^{(m)})+\\varepsilon_{r} k_{0}^{2}\n", "\\mathbf{E}^{(m)}_s \\cdot \\bar{\\mathbf{v}}^{(m)}\n", "+k_{0}^{2}\\left(\\varepsilon_{r}\n", "-\\varepsilon_b\\right)\\mathbf{E}^{(m)}_b \\cdot \\bar{\\mathbf{v}}^{(m)}\\\\\n", "&+\\left(\\boldsymbol{\\mu}^{-1}_{pml} \\nabla \\times \\mathbf{E}^{(m)}_s\n", "\\right)\\cdot \\nabla \\times \\bar{\\mathbf{v}}^{(m)}-k_{0}^{2}\n", "\\left(\\boldsymbol{\\varepsilon}_{pml} \\mathbf{E}^{(m)}_s \\right)\\cdot\n", "\\bar{\\mathbf{v}}^{(m)}~ \\rho d\\rho dz =0\n", "\\end{align}\n", "$$\n", "\n", "In the end, we have multiple weak forms corresponding to the\n", "different cylindrical harmonics, where the integration is performed\n", "over a 2D domain.\n", "\n", "Let's now implement this problem in DOLFINx.\n", "As a first step we can define the function for the $\\nabla\\times$\n", "operator in cylindrical coordinates:" ] }, { "cell_type": "code", "execution_count": null, "id": "d88f56f7", "metadata": {}, "outputs": [], "source": [ "def curl_axis(a, m: int, rho):\n", "\n", " curl_r = -a[2].dx(1) - 1j * m / rho * a[1]\n", " curl_z = a[2] / rho + a[2].dx(0) + 1j * m / rho * a[0]\n", " curl_p = a[0].dx(1) - a[1].dx(0)\n", "\n", " return ufl.as_vector((curl_r, curl_z, curl_p))" ] }, { "cell_type": "markdown", "id": "e70bee2a", "metadata": {}, "source": [ "Then we need to define the analytical formula for the background field.\n", "For our purposes, we can consider the wavevector and the electric field\n", "lying in the same plane of our 2D domain, while the magnetic field is\n", "transverse to such domain. For this reason, we will refer to this polarization\n", "as TMz polarization.\n", "\n", "For a TMz polarization, the cylindrical harmonics $\\mathbf{E}^{(m)}_b$\n", "of the background field can be written in this way\n", "([Wait 1955](https://doi.org/10.1139/p55-024)):\n", "\n", "$$\n", "\\begin{align}\n", "\\mathbf{E}^{(m)}_b = &\\hat{\\rho} \\left(E_{0} \\cos \\theta\n", "e^{i k z \\cos \\theta} i^{-m+1} J_{m}^{\\prime}\\left(k_{0} \\rho \\sin\n", "\\theta\\right)\\right)\\\\\n", "+&\\hat{z} \\left(E_{0} \\sin \\theta e^{i k z \\cos \\theta}i^{-m} J_{m}\n", "\\left(k \\rho \\sin \\theta\\right)\\right)\\\\\n", "+&\\hat{\\phi} \\left(\\frac{E_{0} \\cos \\theta}{k \\rho \\sin \\theta}\n", "e^{i k z \\cos \\theta} i^{-m} J_{m}\\left(k \\rho \\sin \\theta\\right)\\right)\n", "\\end{align}\n", "$$\n", "\n", "with $k = 2\\pi n_b/\\lambda = k_0n_b$ being the wavevector, $\\theta$ being\n", "the angle between $\\mathbf{E}_b$ and $\\hat{\\rho}$, and $J_m$\n", "representing the $m$-th order Bessel function of first kind and\n", "$J_{m}^{\\prime}$ its derivative.\n", "In DOLFINx, we can implement these functions in this way:" ] }, { "cell_type": "code", "execution_count": null, "id": "1416486d", "metadata": {}, "outputs": [], "source": [ "def background_field_rz(theta: float, n_bkg: float, k0: float, m: int, x):\n", "\n", " k = k0 * n_bkg\n", "\n", " a_r = (np.cos(theta) * np.exp(1j * k * x[1] * np.cos(theta))\n", " * (1j)**(-m + 1) * jvp(m, k * x[0] * np.sin(theta), 1))\n", "\n", " a_z = (np.sin(theta) * np.exp(1j * k * x[1] * np.cos(theta))\n", " * (1j)**-m * jv(m, k * x[0] * np.sin(theta)))\n", "\n", " return (a_r, a_z)\n", "\n", "\n", "def background_field_p(theta: float, n_bkg: float, k0: float, m: int, x):\n", "\n", " k = k0 * n_bkg\n", "\n", " a_p = (np.cos(theta) / (k * x[0] * np.sin(theta))\n", " * np.exp(1j * k * x[1] * np.cos(theta)) * m\n", " * (1j)**(-m) * jv(m, k * x[0] * np.sin(theta)))\n", "\n", " return a_p" ] }, { "cell_type": "markdown", "id": "7deb3ee4", "metadata": {}, "source": [ "For PML, we can introduce them in our original domain as a spherical shell.\n", "We can then implement a complex coordinate transformation\n", "of this form in this spherical shell:\n", "\n", "$$\n", "\\begin{align}\n", "& \\rho^{\\prime} = \\rho\\left[1 +j \\alpha/k_0 \\left(\\frac{r\n", "- r_{dom}}{r~r_{pml}}\\right)\\right] \\\\\n", "& z^{\\prime} = z\\left[1 +j \\alpha/k_0 \\left(\\frac{r\n", "- r_{dom}}{r~r_{pml}}\\right)\\right] \\\\\n", "& \\phi^{\\prime} = \\phi \\\\\n", "\\end{align}\n", "$$\n", "\n", "with $\\alpha$ being a parameter tuning the absorption inside the PML,\n", "and $r = \\sqrt{\\rho^2 + z^2}$.\n", "This coordinate transformation has the following jacobian:\n", "\n", "$$\n", "\\mathbf{J}=\\mathbf{A}^{-1}= \\nabla\\boldsymbol{\\rho}^\n", "\\prime(\\boldsymbol{\\rho}) =\n", "\\left[\\begin{array}{ccc}\n", "\\frac{\\partial \\rho^{\\prime}}{\\partial \\rho} &\n", "\\frac{\\partial z^{\\prime}}{\\partial \\rho} &\n", "0 \\\\\n", "\\frac{\\partial \\rho^{\\prime}}{\\partial z} &\n", "\\frac{\\partial z^{\\prime}}{\\partial z} &\n", "0 \\\\\n", "0 &\n", "0 &\n", "\\frac{\\rho^\\prime}{\\rho}\\frac{\\partial \\phi^{\\prime}}{\\partial \\phi}\n", "\\end{array}\\right]=\\left[\\begin{array}{ccc}\n", "J_{11} & J_{12} & 0 \\\\\n", "J_{21} & J_{22} & 0 \\\\\n", "0 & 0 & J_{33}\n", "\\end{array}\\right]\n", "$$\n", "\n", "which we can use to calculate ${\\boldsymbol{\\varepsilon}_{pml}}$ and\n", "${\\boldsymbol{\\mu}_{pml}}$:\n", "\n", "$$\n", "\\begin{align}\n", "& {\\boldsymbol{\\varepsilon}_{pml}} =\n", "A^{-1} \\mathbf{A} {\\boldsymbol{\\varepsilon}_b}\\mathbf{A}^{T}\\\\\n", "& {\\boldsymbol{\\mu}_{pml}} =\n", "A^{-1} \\mathbf{A} {\\boldsymbol{\\mu}_b}\\mathbf{A}^{T}\n", "\\end{align}\n", "$$\n", "\n", "For doing these calculations, we define the\n", "`pml_coordinate` and `create_mu_eps` functions:" ] }, { "cell_type": "code", "execution_count": null, "id": "65485b49", "metadata": {}, "outputs": [], "source": [ "def pml_coordinate(\n", " x, r, alpha: float, k0: float, radius_dom: float, radius_pml: float):\n", "\n", " return (x + 1j * alpha / k0 * x * (r - radius_dom) / (radius_pml * r))\n", "\n", "\n", "def create_eps_mu(pml, rho, eps_bkg, mu_bkg):\n", "\n", " J = ufl.grad(pml)\n", "\n", " # Transform the 2x2 Jacobian into a 3x3 matrix.\n", " J = ufl.as_matrix(((J[0, 0], J[0, 1], 0),\n", " (J[1, 0], J[1, 1], 0),\n", " (0, 0, pml[0] / rho)))\n", "\n", " A = ufl.inv(J)\n", " eps_pml = ufl.det(J) * A * eps_bkg * ufl.transpose(A)\n", " mu_pml = ufl.det(J) * A * mu_bkg * ufl.transpose(A)\n", " return eps_pml, mu_pml" ] }, { "cell_type": "markdown", "id": "9ef4e819", "metadata": {}, "source": [ "We can now define some constants and geometrical parameters,\n", "and then we can generate the mesh with Gmsh, by using the\n", "function `generate_mesh_sphere_axis` in `mesh_sphere_axis.py`:" ] }, { "cell_type": "code", "execution_count": null, "id": "45e4413f", "metadata": {}, "outputs": [], "source": [ "# Constants\n", "epsilon_0 = 8.8541878128 * 10**-12 # Vacuum permittivity\n", "mu_0 = 4 * np.pi * 10**-7 # Vacuum permeability\n", "Z0 = np.sqrt(mu_0 / epsilon_0) # Vacuum impedance\n", "I0 = 0.5 / Z0 # Intensity of electromagnetic field\n", "\n", "# Radius of the sphere\n", "radius_sph = 0.025\n", "\n", "# Radius of the domain\n", "radius_dom = 1\n", "\n", "# Radius of the boundary where scattering efficiency\n", "# is calculated\n", "radius_scatt = 0.4 * radius_dom\n", "\n", "# Radius of the PML shell\n", "radius_pml = 0.25\n", "\n", "# Mesh sizes\n", "mesh_factor = 1\n", "in_sph_size = mesh_factor * 2.0e-3\n", "on_sph_size = mesh_factor * 2.0e-3\n", "scatt_size = mesh_factor * 60.0e-3\n", "pml_size = mesh_factor * 40.0e-3\n", "\n", "# Tags for the subdomains\n", "au_tag = 1\n", "bkg_tag = 2\n", "pml_tag = 3\n", "scatt_tag = 4\n", "\n", "# Mesh generation\n", "model = generate_mesh_sphere_axis(\n", " radius_sph, radius_scatt, radius_dom, radius_pml,\n", " in_sph_size, on_sph_size, scatt_size, pml_size,\n", " au_tag, bkg_tag, pml_tag, scatt_tag)\n", "\n", "domain, cell_tags, facet_tags = model_to_mesh(\n", " model, MPI.COMM_WORLD, 0, gdim=2)\n", "\n", "gmsh.finalize()\n", "MPI.COMM_WORLD.barrier()" ] }, { "cell_type": "markdown", "id": "3cb3c35b", "metadata": {}, "source": [ "Let's have a visual check of the mesh and of the subdomains\n", "by plotting them with PyVista:" ] }, { "cell_type": "code", "execution_count": null, "id": "84e3c01b", "metadata": {}, "outputs": [], "source": [ "if have_pyvista:\n", " topology, cell_types, geometry = plot.create_vtk_mesh(domain, 2)\n", " grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)\n", " pyvista.set_jupyter_backend(\"pythreejs\")\n", " plotter = pyvista.Plotter()\n", " num_local_cells = domain.topology.index_map(domain.topology.dim).size_local\n", " grid.cell_data[\"Marker\"] = \\\n", " cell_tags.values[cell_tags.indices < num_local_cells]\n", " grid.set_active_scalars(\"Marker\")\n", " plotter.add_mesh(grid, show_edges=True)\n", " plotter.view_xy()\n", " if not pyvista.OFF_SCREEN:\n", " plotter.show()\n", " else:\n", " pyvista.start_xvfb()\n", " figure = plotter.screenshot(\"sphere_axis_mesh.png\",\n", " window_size=[1000, 1000])" ] }, { "cell_type": "markdown", "id": "6d96d1b6", "metadata": {}, "source": [ "We can now define our function space. For the $\\hat{\\rho}$ and $\\hat{z}$\n", "components of the electric field, we will use Nedelec elements,\n", "while for the $\\hat{\\phi}$ components we will use Lagrange elements:" ] }, { "cell_type": "code", "execution_count": null, "id": "ec616333", "metadata": {}, "outputs": [], "source": [ "degree = 3\n", "curl_el = ufl.FiniteElement(\"N1curl\", domain.ufl_cell(), degree)\n", "lagr_el = ufl.FiniteElement(\"Lagrange\", domain.ufl_cell(), degree)\n", "V = fem.FunctionSpace(domain, ufl.MixedElement([curl_el, lagr_el]))" ] }, { "cell_type": "markdown", "id": "19628961", "metadata": {}, "source": [ "Let's now define our integration domains:" ] }, { "cell_type": "code", "execution_count": null, "id": "9beb90ea", "metadata": {}, "outputs": [], "source": [ "# Measures for subdomains\n", "dx = ufl.Measure(\"dx\", domain, subdomain_data=cell_tags,\n", " metadata={'quadrature_degree': 60})\n", "\n", "dDom = dx((au_tag, bkg_tag))\n", "dPml = dx(pml_tag)" ] }, { "cell_type": "markdown", "id": "c89d4a7a", "metadata": {}, "source": [ "Let's now define the $\\varepsilon_r$ function:" ] }, { "cell_type": "code", "execution_count": null, "id": "fd7f6c62", "metadata": {}, "outputs": [], "source": [ "n_bkg = 1 # Background refractive index\n", "eps_bkg = n_bkg**2 # Background relative permittivity\n", "eps_au = -1.0782 + 1j * 5.8089\n", "\n", "D = fem.FunctionSpace(domain, (\"DG\", 0))\n", "eps = fem.Function(D)\n", "au_cells = cell_tags.find(au_tag)\n", "bkg_cells = cell_tags.find(bkg_tag)\n", "eps.x.array[au_cells] = np.full_like(\n", " au_cells, eps_au, dtype=np.complex128)\n", "eps.x.array[bkg_cells] = np.full_like(bkg_cells, eps_bkg, dtype=np.complex128)\n", "eps.x.scatter_forward()" ] }, { "cell_type": "markdown", "id": "43c2c876", "metadata": {}, "source": [ "We can now define the characteristic parameters of our background field:" ] }, { "cell_type": "code", "execution_count": null, "id": "fdc0b255", "metadata": {}, "outputs": [], "source": [ "wl0 = 0.4 # Wavelength of the background field\n", "k0 = 2 * np.pi / wl0 # Wavevector of the background field\n", "theta = np.pi / 4 # Angle of incidence of the background field\n", "m_list = [0, 1] # list of harmonics" ] }, { "cell_type": "markdown", "id": "027b2528", "metadata": {}, "source": [ "with `m_list` being a list containing the harmonic\n", "numbers we want to solve the problem for. For\n", "subwavelength structure (as in our case), we can limit\n", "the calculation to few harmonic numbers, i.e.\n", "$m = -1, 0, 1$. Besides, due to the symmetry\n", "of Bessel functions, the solutions for $m = \\pm 1$ are\n", "the same, and therefore we can just consider positive\n", "harmonic numbers." ] }, { "cell_type": "markdown", "id": "17684b4b", "metadata": {}, "source": [ "We can now define `eps_pml` and `mu_pml`:" ] }, { "cell_type": "code", "execution_count": null, "id": "5785941b", "metadata": {}, "outputs": [], "source": [ "rho, z = ufl.SpatialCoordinate(domain)\n", "alpha = 5\n", "r = ufl.sqrt(rho**2 + z**2)\n", "\n", "pml_coords = ufl.as_vector((\n", " pml_coordinate(rho, r, alpha, k0, radius_dom, radius_pml),\n", " pml_coordinate(z, r, alpha, k0, radius_dom, radius_pml)))\n", "\n", "eps_pml, mu_pml = create_eps_mu(pml_coords, rho, eps_bkg, 1)" ] }, { "cell_type": "markdown", "id": "c655e7ae", "metadata": {}, "source": [ "We can now define other objects that will be used inside our\n", "solver loop:" ] }, { "cell_type": "code", "execution_count": null, "id": "1082c577", "metadata": {}, "outputs": [], "source": [ "# Function spaces for saving the solution\n", "V_dg = fem.VectorFunctionSpace(domain, (\"DG\", degree))\n", "V_lagr = fem.FunctionSpace(domain, (\"Lagrange\", degree))\n", "\n", "# Function for saving the rho and z component of the electric field\n", "Esh_rz_m_dg = fem.Function(V_dg)\n", "\n", "# Total field\n", "Eh_m = fem.Function(V)\n", "Esh = fem.Function(V)\n", "\n", "n = ufl.FacetNormal(domain)\n", "n_3d = ufl.as_vector((n[0], n[1], 0))\n", "\n", "# Geometrical cross section of the sphere, for efficiency calculation\n", "gcs = np.pi * radius_sph**2\n", "\n", "# Marker functions for the scattering efficiency integral\n", "marker = fem.Function(D)\n", "scatt_facets = facet_tags.find(scatt_tag)\n", "incident_cells = mesh.compute_incident_entities(domain, scatt_facets,\n", " domain.topology.dim - 1,\n", " domain.topology.dim)\n", "midpoints = mesh.compute_midpoints(\n", " domain, domain.topology.dim, incident_cells)\n", "inner_cells = incident_cells[(midpoints[:, 0]**2\n", " + midpoints[:, 1]**2) < (radius_scatt)**2]\n", "marker.x.array[inner_cells] = 1\n", "\n", "# Define integration domain for the gold sphere\n", "dAu = dx(au_tag)\n", "\n", "# Define integration facet for the scattering efficiency\n", "dS = ufl.Measure(\"dS\", domain, subdomain_data=facet_tags)\n", "\n", "phi = 0" ] }, { "cell_type": "markdown", "id": "49e27a24", "metadata": {}, "source": [ "We can now solve our problem for all the chosen harmonic numbers:" ] }, { "cell_type": "code", "execution_count": null, "id": "1fc3aa95", "metadata": {}, "outputs": [], "source": [ "for m in m_list:\n", "\n", " # Definition of Trial and Test functions\n", " Es_m = ufl.TrialFunction(V)\n", " v_m = ufl.TestFunction(V)\n", "\n", " # Background field\n", " Eb_m = fem.Function(V)\n", " f_rz = partial(background_field_rz, theta, n_bkg, k0, m)\n", " f_p = partial(background_field_p, theta, n_bkg, k0, m)\n", " Eb_m.sub(0).interpolate(f_rz)\n", " Eb_m.sub(1).interpolate(f_p)\n", "\n", " curl_Es_m = curl_axis(Es_m, m, rho)\n", " curl_v_m = curl_axis(v_m, m, rho)\n", "\n", " F = - ufl.inner(curl_Es_m, curl_v_m) * rho * dDom \\\n", " + eps * k0 ** 2 * ufl.inner(Es_m, v_m) * rho * dDom \\\n", " + k0 ** 2 * (eps - eps_bkg) * ufl.inner(Eb_m, v_m) * rho * dDom \\\n", " - ufl.inner(ufl.inv(mu_pml) * curl_Es_m, curl_v_m) * rho * dPml \\\n", " + k0 ** 2 * ufl.inner(eps_pml * Es_m, v_m) * rho * dPml\n", "\n", " a, L = ufl.lhs(F), ufl.rhs(F)\n", "\n", " problem = fem.petsc.LinearProblem(a, L, bcs=[], petsc_options={\n", " \"ksp_type\": \"preonly\", \"pc_type\": \"lu\"})\n", " Esh_m = problem.solve()\n", "\n", " Esh_rz_m, Esh_p_m = Esh_m.split()\n", "\n", " # Interpolate over rho and z components over DG space\n", " Esh_rz_m_dg.interpolate(Esh_rz_m)\n", "\n", " # Save solutions\n", " with VTXWriter(domain.comm, f\"sols/Es_rz_{m}.bp\", Esh_rz_m_dg) as f:\n", " f.write(0.0)\n", " with VTXWriter(domain.comm, f\"sols/Es_p_{m}.bp\", Esh_p_m) as f:\n", " f.write(0.0)\n", "\n", " # Define scattered magnetic field\n", " Hsh_m = 1j * curl_axis(Esh_m, m, rho) / (Z0 * k0 * n_bkg)\n", "\n", " # Total electric field\n", " Eh_m.x.array[:] = Eb_m.x.array[:] + Esh_m.x.array[:]\n", "\n", " if m == 0:\n", "\n", " Esh.x.array[:] = Esh_m.x.array[:] * np.exp(- 1j * m * phi)\n", "\n", " elif m == m_list[0]:\n", "\n", " Esh.x.array[:] = 2 * Esh_m.x.array[:] * np.exp(- 1j * m * phi)\n", "\n", " else:\n", "\n", " Esh.x.array[:] += 2 * Esh_m.x.array[:] * np.exp(- 1j * m * phi)\n", "\n", " # Efficiencies calculation\n", "\n", " if m == 0: # initialize and do not add 2 factor\n", " P = np.pi * ufl.inner(ufl.cross(Esh_m, ufl.conj(Hsh_m)), n_3d) * marker\n", " Q = np.pi * eps_au.imag * k0 * (ufl.inner(Eh_m, Eh_m)) / Z0 / n_bkg\n", "\n", " q_abs_fenics_proc = (fem.assemble_scalar(\n", " fem.form(Q * rho * dAu)) / gcs / I0).real\n", " q_sca_fenics_proc = (fem.assemble_scalar(\n", " fem.form((P('+') + P('-')) * rho * dS(scatt_tag))) / gcs / I0).real\n", " q_abs_fenics = domain.comm.allreduce(q_abs_fenics_proc, op=MPI.SUM)\n", " q_sca_fenics = domain.comm.allreduce(q_sca_fenics_proc, op=MPI.SUM)\n", "\n", " elif m == m_list[0]: # initialize and add 2 factor\n", " P = 2 * np.pi * ufl.inner(ufl.cross(Esh_m,\n", " ufl.conj(Hsh_m)), n_3d) * marker\n", " Q = 2 * np.pi * eps_au.imag * k0 * (ufl.inner(Eh_m, Eh_m)) / Z0 / n_bkg\n", "\n", " q_abs_fenics_proc = (fem.assemble_scalar(\n", " fem.form(Q * rho * dAu)) / gcs / I0).real\n", " q_sca_fenics_proc = (fem.assemble_scalar(\n", " fem.form((P('+') + P('-')) * rho * dS(scatt_tag))) / gcs / I0).real\n", " q_abs_fenics = domain.comm.allreduce(q_abs_fenics_proc, op=MPI.SUM)\n", " q_sca_fenics = domain.comm.allreduce(q_sca_fenics_proc, op=MPI.SUM)\n", "\n", " else: # do not initialize and add 2 factor\n", " P = 2 * np.pi * ufl.inner(ufl.cross(Esh_m,\n", " ufl.conj(Hsh_m)), n_3d) * marker\n", " Q = 2 * np.pi * eps_au.imag * k0 * (ufl.inner(Eh_m, Eh_m)) / Z0 / n_bkg\n", "\n", " q_abs_fenics_proc = (fem.assemble_scalar(\n", " fem.form(Q * rho * dAu)) / gcs / I0).real\n", " q_sca_fenics_proc = (fem.assemble_scalar(\n", " fem.form((P('+') + P('-')) * rho * dS(scatt_tag))) / gcs / I0).real\n", " q_abs_fenics += domain.comm.allreduce(q_abs_fenics_proc, op=MPI.SUM)\n", " q_sca_fenics += domain.comm.allreduce(q_sca_fenics_proc, op=MPI.SUM)\n", "\n", "q_ext_fenics = q_abs_fenics + q_sca_fenics" ] }, { "cell_type": "markdown", "id": "f77b41c9", "metadata": {}, "source": [ "Let's compare the analytical and numerical efficiencies, and let's print\n", "the results:" ] }, { "cell_type": "code", "execution_count": null, "id": "4c76319f", "metadata": {}, "outputs": [], "source": [ "q_abs_analyt = 0.9622728008329892\n", "q_sca_analyt = 0.07770397394691526\n", "q_ext_analyt = q_abs_analyt + q_sca_analyt" ] }, { "cell_type": "code", "execution_count": null, "id": "4a97c5d7", "metadata": {}, "outputs": [], "source": [ "# Error calculation\n", "err_abs = np.abs(q_abs_analyt - q_abs_fenics) / q_abs_analyt\n", "err_sca = np.abs(q_sca_analyt - q_sca_fenics) / q_sca_analyt\n", "err_ext = np.abs(q_ext_analyt - q_ext_fenics) / q_ext_analyt\n", "\n", "if MPI.COMM_WORLD.rank == 0:\n", "\n", " print()\n", " print(f\"The analytical absorption efficiency is {q_abs_analyt}\")\n", " print(f\"The numerical absorption efficiency is {q_abs_fenics}\")\n", " print(f\"The error is {err_abs*100}%\")\n", " print()\n", " print(f\"The analytical scattering efficiency is {q_sca_analyt}\")\n", " print(f\"The numerical scattering efficiency is {q_sca_fenics}\")\n", " print(f\"The error is {err_sca*100}%\")\n", " print()\n", " print(f\"The analytical extinction efficiency is {q_ext_analyt}\")\n", " print(f\"The numerical extinction efficiency is {q_ext_fenics}\")\n", " print(f\"The error is {err_ext*100}%\")\n", "\n", " assert err_abs < 0.01\n", " assert err_sca < 0.01\n", " assert err_ext < 0.01\n", "\n", "\n", "Esh_rz, Esh_p = Esh.split()\n", "\n", "Esh_rz_dg = fem.Function(V_dg)\n", "Esh_r_dg = fem.Function(V_dg)\n", "\n", "# Interpolate over rho and z components over DG space\n", "Esh_rz_dg.interpolate(Esh_rz)\n", "\n", "with VTXWriter(domain.comm, \"sols/Es_rz.bp\", Esh_rz_dg) as f:\n", " f.write(0.0)\n", "with VTXWriter(domain.comm, \"sols/Es_p.bp\", Esh_p) as f:\n", " f.write(0.0)\n", "\n", "if have_pyvista:\n", " V_cells, V_types, V_x = plot.create_vtk_mesh(V_dg)\n", " V_grid = pyvista.UnstructuredGrid(V_cells, V_types, V_x)\n", " Esh_r_values = np.zeros((V_x.shape[0], 3), dtype=np.float64)\n", " Esh_r_values[:, :domain.topology.dim] = \\\n", " Esh_rz_dg.x.array.reshape(V_x.shape[0], domain.topology.dim).real\n", "\n", " V_grid.point_data[\"u\"] = Esh_r_values\n", "\n", " pyvista.set_jupyter_backend(\"pythreejs\")\n", " plotter = pyvista.Plotter()\n", "\n", " plotter.add_text(\"magnitude\", font_size=12, color=\"black\")\n", " plotter.add_mesh(V_grid.copy(), show_edges=False)\n", " plotter.view_xy()\n", " plotter.link_views()\n", "\n", " if not pyvista.OFF_SCREEN:\n", " plotter.show()\n", " else:\n", " pyvista.start_xvfb()\n", " plotter.screenshot(\"Esh_r.png\", window_size=[800, 800])" ] } ], "metadata": { "jupytext": { "formats": "ipynb,py:light" }, "kernelspec": { "display_name": "Python 3 (DOLFINx complex)", "language": "python", "name": "python3-complex" }, "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.10.4" } }, "nbformat": 4, "nbformat_minor": 5 }