complete sheet4 exercise 4.7

campospinto · Mar 11, 2025 · Mar 11, 2025 · Mar 24, 2025 · Mar 24, 2025 · Mar 24, 2025
commit f9537c01b7f07ca100a803f2429a28ec06b3fe59
diff --git a/exercises/hodge_laplace_annulus.ipynb b/exercises/hodge_laplace_annulus.ipynb
@@ -1,5 +1,71 @@
 {
  "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "**Note: This notebook must be changed by the student. As it is now, it does not approximate the solution to the Poisson problem!**"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# 2D Solver for the Magneto-Static Problem on an Annulus\n",
+    "\n",
+    "In this exercise we write a solver for the 2D magneto-static problem on an annulus (r_min=.3, r_max=1.) $\\Omega$.\n",
+    "\n",
+    "\\begin{align*}\n",
+    "    &\\text{Given }\\boldsymbol{J}\\in H_0(curl;\\Omega)\\cap(\\mathfrak{H}^1)^{\\perp}\\text{ s.th. }div\\ \\boldsymbol{J}=0, \\\\\n",
+    "    &\\text{find }B\\in H(\\boldsymbol{curl};\\Omega)\\text{ satisfying} \\\\\n",
+    "    &\\boldsymbol{curl}\\ B = \\boldsymbol{J} \\quad \\text{ in }\\Omega,\n",
+    "\\end{align*}\n",
+    "\n",
+    "for the specific choice\n",
+    "\n",
+    "\\begin{align*}\n",
+    "    r(x, y) &= \\sqrt{x^2 + y^2}, \\\\\n",
+    "    B(x, y) &= \\frac{-x}{r^3}, \\\\\n",
+    "    \\boldsymbol{J}(x, y) &= \\boldsymbol{curl}\\ B.\n",
+    "\\end{align*}\n",
+    "\n",
+    "Note that $B$ is chosen to be the $curl$ of $\\boldsymbol{A}$ such that $\\boldsymbol{A}\\in H_0(curl;\\Omega)$ and $div\\ \\boldsymbol{A} = 0$ for\n",
+    "\\begin{align}\n",
+    "    \\boldsymbol{A}(x, y) = \\frac{1}{r^3}\\left(\\begin{matrix} xy \\\\ y^2 \\end{matrix}\\right).\n",
+    "\\end{align}                \n",
+    "\n",
+    "This problem is ill-posed. We have shown, that the following auxiliary problem is well-posed, and that its solution $\\boldsymbol{A}$ satisfies $\\boldsymbol{curl}\\ curl\\ \\boldsymbol{A} = \\boldsymbol{J}$, hence $B = curl\\ \\boldsymbol{A}$ is a solution to the original problem.\n",
+    "\n",
+    "\\begin{align*}\n",
+    "    &\\text{Given }\\boldsymbol{J}\\in H_0(curl;\\Omega)\\cap(\\mathfrak{H}^1)^{\\perp}\\text{ s.th. }div\\ \\boldsymbol{J}=0, \\\\\n",
+    "    &\\text{find }\\boldsymbol{A}\\in D(\\mathcal{L}^1)\\cap(\\mathfrak{H}^1)^{\\perp}\\text{ satisfying} \\\\\n",
+    "    &\\mathcal{L}^1(\\boldsymbol{A}) = \\boldsymbol{J} \\quad \\text{ in }\\Omega.\n",
+    "\\end{align*}\n",
+    "\n",
+    "## The Discrete Variational Formulation\n",
+    "\n",
+    "The corresponding discrete variational formulation reads\n",
+    "\n",
+    "\\begin{align*}\n",
+    "    &\\quad\\text{For }\\boldsymbol{J}\\in V_h^1\\cap(\\mathfrak{H}_h^1)^{\\perp}\\text{ s.th. }\\widetilde{div}_h\\ \\boldsymbol{J}=0, \\\\\n",
+    "    &\\quad\\text{find } (\\sigma,\\boldsymbol{A},\\boldsymbol{p})\\in V_h^0\\times V_h^1\\times\\mathfrak{H}_h^1 \\text{ satisfying} \\\\\n",
+    "    (\\textbf{HL1}_h)\n",
+    "    &\\begin{cases}\n",
+    "        \\quad\\quad\\quad\\ \\ \\langle\\sigma,\\tau\\rangle - \\quad\\ \\ \\langle\\boldsymbol{A},\\boldsymbol{grad}\\ \\tau\\rangle &= 0 \\ \\quad\\quad\\forall\\tau\\in V_h^0, \\\\\n",
+    "        -\\langle\\boldsymbol{grad}\\ \\sigma,\\boldsymbol{v}\\rangle + \\langle curl\\ \\boldsymbol{A},curl\\ \\boldsymbol{v}\\rangle + \\langle\\boldsymbol{p},\\boldsymbol{v}\\rangle &= \\langle\\boldsymbol{J},\\boldsymbol{v}\\rangle \\ \\forall\\boldsymbol{v}\\in V_h^1, \\\\\n",
+    "        \\quad\\quad\\quad\\quad\\quad\\quad\\quad\\quad\\quad\\quad\\quad\\langle\\boldsymbol{A},\\boldsymbol{q}\\rangle &= 0 \\ \\quad\\quad\\forall \\boldsymbol{q}\\in\\mathfrak{H}_h^1.\n",
+    "    \\end{cases}\n",
+    "\\end{align*}\n",
+    "\n",
+    "The discrete solution $\\boldsymbol{A}$ will then satisfy $\\widetilde{\\boldsymbol{curl}}_h\\ curl\\ \\boldsymbol{A} = \\boldsymbol{J}$, hence $B = curl\\ \\boldsymbol{A}$ will satisfy $\\widetilde{\\boldsymbol{curl}}_h B = \\boldsymbol{J}$.\n",
+    "\n",
+    "## Method of manufactured solution\n",
+    "\n",
+    "To verify that our code works, we employ the method of manufactured solutions. Given the vector field $\\boldsymbol{A}_{ex}$, see (1), we compute the r.h.s. via $\\langle\\boldsymbol{J},\\ \\boldsymbol{v}\\rangle = \\langle\\widetilde{\\boldsymbol{curl}}_h\\ curl\\ \\boldsymbol{A}_{ex},\\ \\boldsymbol{v}\\rangle = \\langle curl\\ \\boldsymbol{A}_{ex},\\ curl\\ \\boldsymbol{v}\\rangle$.\n",
+    "\n",
+    "We then report the $L^2$ errors $\\| \\boldsymbol{A}_h - \\boldsymbol{A}_{ex} \\|$ and $\\| B_h - B_{ex} \\|$."
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -46,8 +112,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "from sympde.calculus            import dot\n",
-    "from sympde.expr                import LinearForm, BilinearForm, integral\n",
+    "from sympde.calculus            import dot, curl\n",
+    "from sympde.expr                import LinearForm, BilinearForm, integral, Norm\n",
     "from sympde.topology            import elements_of, Derham\n",
     "\n",
     "derham = Derham(domain, sequence=['h1', 'hcurl', 'l2'])\n",
@@ -168,7 +234,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## Boundary Conditions and inverse Matrices\n",
+    "## Boundary Conditions and Inverse Matrices\n",
     "\n",
     "Let us denote the interior coefficients by a subscript $I$ and the exterior coefficients by a subscript $\\Gamma$. Suppose we want to solve the linear system\n",
     "$$\n",
@@ -208,7 +274,7 @@
     "from psydac.linalg.solvers import inverse\n",
     "\n",
     "M0_0    = P_H1 @ M0 @ P_H1 + P_H1_Gamma\n",
-    "M0_inv  = inverse(M0_0, 'cg', maxiter=1000, tol=1e-9)"
+    "M0_inv  = inverse(M0_0, 'cg', maxiter=1000, tol=1e-12)"
    ]
   },
   {
@@ -220,8 +286,8 @@
     "Psydac does not offer a direct way to access eigenvalues and -vectors of matrices.\n",
     "However, we can transform our operators to scipy.sparse format and use scipy's `eigsh` method.\n",
     "\n",
-    "We want our solution vector field $A$ to be an element of $(\\mathfrak{h}_h^1)^{\\perp}$, \n",
-    "and it holds $\\ker\\mathcal{L}_h^1 = \\mathfrak{h}_h^1$.\n",
+    "We want our solution vector field $A$ to be an element of $(\\mathfrak{H}_h^1)^{\\perp}$, \n",
+    "and it holds $\\ker\\mathcal{L}_h^1 = \\mathfrak{H}_h^1$.\n",
     "\n",
     "Because the annulus has one whole, and hence the space of harmonic forms is of dimension 1, we want to compute the first eigenvector (corresponding to the eigenvalue $0$) of the operator\n",
     "$$\n",
@@ -271,9 +337,9 @@
     "inv_M0_sp.eliminate_zeros()\n",
     "\n",
     "# For this complex operator, the projection method consists of using the modified gradient\n",
-    "L1_sp               = Ct_sp @ M2_sp @ C_sp + M1_sp @ G_sp_0 @ inv_M0_sp @ Gt_sp_0 @ M1_sp\n",
+    "#L1_sp               =                                                                          # <- todo \n",
     "# ... and the projection matrices as usual\n",
-    "L1_sp_bc            = P_Hcurl_sp @ L1_sp @ P_Hcurl_sp + P_Hcurl_Gamma_sp"
+    "#L1_sp_bc            = P_Hcurl_sp @ L1_sp @ P_Hcurl_sp + P_Hcurl_Gamma_sp"
    ]
   },
   {
@@ -292,14 +358,14 @@
     "from utils                  import get_eigenvalues\n",
     "\n",
     "dim_harmonic_space          = 1\n",
-    "eigenvalues, eigenvectors   = get_eigenvalues(dim_harmonic_space, 1e-6, L1_sp_bc, M1_sp)\n",
+    "#eigenvalues, eigenvectors   = get_eigenvalues(dim_harmonic_space, 1e-6, L1_sp_bc, M1_sp)\n",
     "\n",
     "# A basis of the vector space of harmonic forms (hf)\n",
-    "hf                  = eigenvectors[:,0]\n",
+    "#hf                  = eigenvectors[:,0]\n",
     "# Obtain the eigenvector as a \"row csc_matrix (1 x n)\"\n",
-    "hft_sp              = csc_matrix(hf)\n",
+    "#hft_sp              = csc_matrix(hf)\n",
     "# ... and as a \"column csc_matrix (n x 1)\"\n",
-    "hf_sp               = hft_sp.transpose()"
+    "#hf_sp               = hft_sp.transpose()"
    ]
   },
   {
@@ -317,9 +383,9 @@
    "outputs": [],
    "source": [
     "# System matrix without taking BCs into account\n",
-    "A_sp = bmat([[M0_sp,            Gt_sp @ M1_sp,          None],\n",
-    "             [M1_sp @ G_sp,     Ct_sp @ M2_sp @ C_sp,   M1_sp @ hf_sp],\n",
-    "             [None,             hft_sp @ M1_sp,         None]])\n",
+    "#A_sp = bmat([[x,            x,          None],                                 # <- todo\n",
+    "#             [x,     x,   M1_sp @ x],                                          # <- todo\n",
+    "#             [None,             x,         None]])                             # <- todo\n",
     "\n",
     "ZERO    = csc_matrix(np.array([0]))\n",
     "ONE     = csc_matrix(np.array([1]))\n",
@@ -333,29 +399,7 @@
     "                [None,          None,               ZERO]])\n",
     "\n",
     "# System matrix taking BCs into account\n",
-    "A_sp_bc = P @ A_sp @ P + P_Gamma"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "# Method of manufactured solution\n",
-    "\n",
-    "In order to test our code, we employ the method of manufactured solution.\n",
-    "We start from the true $\\boldsymbol{A}\\in H_0(curl;\\Omega),\\ div\\ \\boldsymbol{A} = 0$\n",
-    "$$\n",
-    "\\boldsymbol{A}(x, y) = \\frac{1}{(x^2 + y^2)^{3/2}}\\left(\\begin{matrix} xy \\\\ y^2 \\end{matrix}\\right)\n",
-    "$$\n",
-    "Because the r.h.s. of the discrete variational formulation reads\n",
-    "$$\n",
-    "r.h.s. = \\left(\\begin{matrix} \\mathbb{0} \\\\ (\\boldsymbol{v},\\ \\boldsymbol{J}) \\\\ \\mathbb{0} \\end{matrix}\\right)\n",
-    "$$\n",
-    "and because $\\boldsymbol{J} = \\boldsymbol{curl}\\ curl\\ \\boldsymbol{A}$, the r.h.s. of the linear system becomes\n",
-    "$$\n",
-    "r.h.s. = \\left(\\begin{matrix} \\mathbb{0} \\\\ \\mathbb{C}^T\\mathbb{M}_2\\mathbb{C}\\mathbb{A} \\\\ \\mathbb{0} \\end{matrix}\\right)\n",
-    "$$\n",
-    "with $\\mathbb{A}$ being the coefficient vector of discrete $\\boldsymbol{A}$. If our code is correct, we will obtain approximately $\\mathbb{A}$ as second component of the solution vector."
+    "#A_sp_bc = P @ A_sp @ P + P_Gamma"
    ]
   },
   {
@@ -364,32 +408,33 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Define exact solution A_ex for method of manufactured solution\n",
-    "r       = lambda x, y : np.sqrt(x**2 + y**2)\n",
-    "A_ex_1  = lambda x, y : (x*y)  / (r(x,y)**3)\n",
-    "A_ex_2  = lambda x, y : (y**2) / (r(x,y)**3)\n",
+    "from sympy import sqrt, Tuple, Matrix\n",
+    "x, y            = domain.coordinates\n",
+    "r_sympy         = sqrt(x**2 + y**2)\n",
+    "A_ex_1_sympy    = (x*y)  / (r_sympy**3)\n",
+    "A_ex_2_sympy    = (y**2) / (r_sympy**3)\n",
+    "A_ex_sympy      = Tuple(A_ex_1_sympy, A_ex_2_sympy)\n",
     "\n",
-    "A_ex            = (A_ex_1, A_ex_2)\n",
-    "A_ex_FemField   = P1(A_ex)\n",
-    "A_ex_coeffs     = A_ex_FemField.coeffs\n",
+    "#l               =                                                      # <- todo\n",
+    "#lh              = discretize(l, domain_h, V1h)\n",
     "\n",
     "# First component of rhs\n",
     "rhs1_nparray    = np.zeros(V0h.nbasis)\n",
     "\n",
     "# Second component of rhs\n",
-    "rhs2_nparray    = (P_Hcurl @ Ct @ M2 @ C @ A_ex_coeffs ).toarray()\n",
+    "#rhs2_nparray    = (P_Hcurl @ lh.assemble()).toarray()\n",
     "\n",
     "# Third component of rhs\n",
     "rhs3_nparray    = np.zeros(dim_harmonic_space)\n",
     "\n",
     "# Full rhs\n",
-    "rhs_nparray     = np.block([rhs1_nparray, rhs2_nparray, rhs3_nparray])\n",
+    "#rhs_nparray     = np.block([rhs1_nparray, rhs2_nparray, rhs3_nparray])\n",
     "\n",
     "# -------------- direct solve with scipy spsolve ---------------\n",
     "import  time\n",
     "from    scipy.sparse.linalg import spsolve\n",
     "t0 = time.time()\n",
-    "sol_nparray = spsolve(A_sp_bc.asformat('csr'), rhs_nparray)\n",
+    "#sol_nparray = spsolve(A_sp_bc.asformat('csr'), rhs_nparray)\n",
     "t1 = time.time()\n",
     "# --------------------------------------------------------------"
    ]
@@ -407,22 +452,84 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "sigmah_nparray  = sol_nparray[:V0h.nbasis]\n",
-    "Ah_nparray      = sol_nparray[V0h.nbasis:V0h.nbasis + V1h.nbasis]\n",
-    "ph_nparray      = sol_nparray[V0h.nbasis+V1h.nbasis:]\n",
+    "# print the result\n",
+    "print( '> Grid              :: [{ne1},{ne2}]'.format( ne1=ncells[0], ne2=ncells[1]) )\n",
+    "print( '> Degree            :: [{p1},{p2}]'  .format( p1=degree[0], p2=degree[1] ) )\n",
+    "print( '> Solution time     :: {:.3g}'.format( t1-t0 ) )\n",
+    "print( '' )\n",
     "\n",
-    "from psydac.linalg.utilities import array_to_psydac\n",
-    "Ah_coeffs       = array_to_psydac(Ah_nparray, V1h.coeff_space)\n",
+    "#sigmah_nparray  = sol_nparray[:V0h.nbasis]\n",
+    "#Ah_nparray      = sol_nparray[V0h.nbasis:V0h.nbasis + V1h.nbasis]\n",
+    "#ph_nparray      = sol_nparray[V0h.nbasis+V1h.nbasis:]\n",
+    "\n",
+    "from psydac.linalg.utilities    import array_to_psydac\n",
+    "from psydac.fem.basic           import FemField\n",
+    "#Ah_coeffs       = array_to_psydac(Ah_nparray, V1h.coeff_space)\n",
+    "\n",
+    "# ----- Ah vs. Aex -----\n",
+    "#Ah_FemField     = FemField(V1h, Ah_coeffs)\n",
+    "error           = Matrix([u1[0] - A_ex_sympy[0], u1[1] - A_ex_sympy[1]])\n",
+    "l2norm          = Norm(error, domain, kind='l2')\n",
+    "l2norm_h        = discretize(l2norm, domain_h, V1h, backend=backend)\n",
+    "#l2_error        = l2norm_h.assemble(u1=Ah_FemField)\n",
+    "\n",
+    "#print( '> || Ah - Aex ||    :: {:.2e}'.format( l2_error ) )\n",
+    "# ----------------------\n",
+    "\n",
+    "# ----- Bh vs. Bex -----\n",
+    "B_ex_sympy      = ( -x ) / ( r_sympy**3 )\n",
+    "#Bh_coeffs       = C @ Ah_coeffs\n",
+    "#Bh_FemField     = FemField(V2h, Bh_coeffs)\n",
+    "error           = u2 - B_ex_sympy\n",
+    "l2norm          = Norm(error, domain, kind='l2')\n",
+    "l2norm_h        = discretize(l2norm, domain_h, V2h, backend=backend)\n",
+    "#l2_error        = l2norm_h.assemble(u2=Bh_FemField)\n",
+    "\n",
+    "#print( '> || Bh - Bex ||    :: {:.2e}'.format( l2_error ) )\n",
+    "# ----------------------"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from sympy                          import lambdify\n",
     "\n",
-    "error           = A_ex_coeffs - Ah_coeffs\n",
-    "l2_error        = np.sqrt( error.dot(M1 @ error) )\n",
+    "from psydac.fem.plotting_utilities  import plot_field_2d\n",
     "\n",
-    "# print the result\n",
-    "print( '> Grid          :: [{ne1},{ne2}]'.format( ne1=ncells[0], ne2=ncells[1]) )\n",
-    "print( '> Degree        :: [{p1},{p2}]'  .format( p1=degree[0], p2=degree[1] ) )\n",
-    "print( '> L2 error      :: {:.2e}'.format( l2_error ) )\n",
-    "print( '' )\n",
-    "print( '> Solution time :: {:.3g}'.format( t1-t0 ) )"
+    "Aex_lambdified  = lambdify((x, y), A_ex_sympy)\n",
+    "#A_error_fun     = lambda x, y : Aex_lambdified(x, y) - Ah_FemField(x, y)\n",
+    "\n",
+    "Bex_lambdified  = lambdify((x, y), B_ex_sympy)\n",
+    "#B_error_fun     = lambda x, y : Bex_lambdified(x, y) - Bh_FemField(x, y)\n",
+    "\n",
+    "#plot_field_2d(Ah_FemField, domain=domain, space_kind='hcurl',   N_vis=250,  title=r'$\\| \\boldsymbol{A}_h \\|$', filename=None, plot_type='vector_field')\n",
+    "#plot_field_2d(Ah_FemField, domain=domain, space_kind='hcurl',   N_vis=50,   title=r'$\\| \\boldsymbol{A}_h \\|$', filename=None)\n",
+    "\n",
+    "#plot_field_2d(Bh_FemField, domain=domain, space_kind='l2',      N_vis=500,  title=r'$\\|        B_h       \\|$', filename=None, surface_plot=True)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Harmonic Field\n",
+    "\n",
+    "We can also plot the harmonic field"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#hf_coeffs   = array_to_psydac(hf, V1h.coeff_space)\n",
+    "\n",
+    "#plot_field_2d(stencil_coeffs=hf_coeffs, Vh=V1h, domain=domain, space_kind='hcurl', N_vis=50,  title=r'$\\| \\boldsymbol{\\mathfrak{H}} \\|$', filename=None)\n",
+    "#plot_field_2d(stencil_coeffs=hf_coeffs, Vh=V1h, domain=domain, space_kind='hcurl', N_vis=250, title=r'$\\| \\boldsymbol{\\mathfrak{H}} \\|$', filename=None, plot_type='vector_field')"
    ]
   }
  ],