diff --git a/.gitignore b/.gitignore index 17ee1da..c25dd40 100644 --- a/.gitignore +++ b/.gitignore @@ -19,3 +19,5 @@ usr .env .iga-python *_autosummary* + +.DS_Store \ No newline at end of file diff --git a/_toc.yml b/_toc.yml index 1f7600b..e1dac24 100644 --- a/_toc.yml +++ b/_toc.yml @@ -66,6 +66,10 @@ parts: - file: chapter3/navier-stokes-steady sections: - file: chapter3/navier-stokes-steady-streamfunction-velocity + - caption: Exercises + chapters: + - file: exercises/harmonic-fem + - file: exercises/harmonic-feec - caption: Advanced subjects chapters: - file: chapter4/subdomains diff --git a/exercises/electro_static_2d.ipynb b/exercises/electro_static_2d.ipynb new file mode 100644 index 0000000..0a7b1d5 --- /dev/null +++ b/exercises/electro_static_2d.ipynb @@ -0,0 +1,285 @@ +{ + "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 Electro-Static Problem on the Unit Square\n", + "\n", + "In this exercise we write a solver for the 2D electro-static problem on the unit square, which can be shown to be equivalent to an $\\mathcal{L}^1$ Hodge-Laplace problem.\n", + "\n", + "\\begin{align*}\n", + " &\\text{Find }\\boldsymbol{E} \\in D(\\mathcal{L}^1) \\\\\n", + " &\\mathcal{L}^1\\boldsymbol{E} = (\\boldsymbol{curl}\\ curl - \\boldsymbol{grad}\\ div)\\boldsymbol{E} = -\\boldsymbol{grad}\\ \\rho \\quad \\text{ in }\\Omega=]0,1[^2,\n", + "\\end{align*}\n", + "for the specific choice\n", + "\\begin{align*}\n", + " \\rho(x, y) = 2\\sin(\\pi x)\\sin(\\pi y)\n", + "\\end{align*}\n", + "\n", + "## The Discrete Variational Formulation\n", + "\n", + "The corresponding discrete variational formulation reads\n", + "\n", + "\\begin{align*}\n", + " &\\text{Find }\\boldsymbol{E} \\in V_h^1 \\subset V^1 = H_0(curl;\\Omega)\\text{ such that } \\\\\n", + " &a(\\boldsymbol{E}, \\boldsymbol{v}) = l(\\boldsymbol{v})\\quad\\forall\\ \\boldsymbol{v}\\in V_h^1\n", + "\\end{align*}\n", + "\n", + "where \n", + "- $a(\\boldsymbol{E},\\boldsymbol{v}) \\coloneqq \\int_{\\Omega} (curl\\ \\boldsymbol{E})(curl\\ \\boldsymbol{v}) + (\\widetilde{div}_h\\boldsymbol{E})(\\widetilde{div}_h\\boldsymbol{v}) ~ d\\Omega$ ,\n", + "- $l(\\boldsymbol{v}) := -\\int_{\\Omega} \\boldsymbol{grad}\\ \\rho\\cdot\\boldsymbol{v} ~ d\\Omega$." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Discrete Model using de Rham objects" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from sympde.calculus import dot\n", + "from sympde.topology import elements_of, Square, Derham\n", + "from sympde.expr.expr import BilinearForm, integral\n", + "\n", + "domain = Square('S', bounds1=(0, 1), bounds2=(0, 1))\n", + "derham = Derham(domain, sequence=['h1', 'hcurl', 'l2'])\n", + "\n", + "V0 = derham.V0\n", + "V1 = derham.V1\n", + "V2 = derham.V2\n", + "\n", + "u0, v0 = elements_of(V0, names='u0, v0')\n", + "u1, v1 = elements_of(V1, names='u1, v1')\n", + "u2, v2 = elements_of(V2, names='u2, v2')\n", + "\n", + "# bilinear forms corresponding to V0, V1 and V2 mass matrices\n", + "m0 = BilinearForm((u0, v0), integral(domain, u0 * v0))\n", + "m1 = BilinearForm((u1, v1), integral(domain, dot(u1, v1)))\n", + "m2 = BilinearForm((u2, v2), integral(domain, u2 * v2))\n", + "\n", + "# callable source term \\rho\n", + "from sympy import pi, sin, cos, lambdify, Tuple, Matrix\n", + "x,y = domain.coordinates\n", + "rho = 2*sin(pi*x)*sin(pi*y)\n", + "rho_lambdified = lambdify((x, y), rho)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from psydac.api.discretization import discretize\n", + "from psydac.api.settings import PSYDAC_BACKEND_GPYCCEL\n", + "\n", + "backend = PSYDAC_BACKEND_GPYCCEL" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ncells = [32, 32] # Bspline cells\n", + "degree = [4, 4] # Bspline degree" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# discretize domain and derham\n", + "# <-\n", + "# <-\n", + "\n", + "# define FEM spaces V0_h, V1_h and V2_h\n", + "#V0_h = derham_h.V0\n", + "#V1_h = derham_h.V1\n", + "#V2_h = derham_h.V2\n", + "\n", + "# Commuting projection operators\n", + "#P0, P1, P2 = derham_h.projectors()\n", + "\n", + "# Exterior Derivative operators (grad and curl)\n", + "# <-\n", + "\n", + "# Mass matrices\n", + "# <-\n", + "# <-\n", + "# <-\n", + "\n", + "# <-\n", + "# <-\n", + "# <-\n", + "\n", + "# Boundary Projectors\n", + "from utils import H1BoundaryProjector2D, HcurlBoundaryProjector2D\n", + "from psydac.linalg.basic import IdentityOperator\n", + "\n", + "#P_H1 = H1BoundaryProjector2D(V0, V0_h.coeff_space)\n", + "#P_Hcurl = HcurlBoundaryProjector2D(V1, V1_h.coeff_space)\n", + "\n", + "#I0 = IdentityOperator(V0_h.coeff_space)\n", + "#I1 = IdentityOperator(V1_h.coeff_space)\n", + "\n", + "#P_H1_Gamma = I0 - P_H1\n", + "#P_Hcurl_Gamma = I1 - P_Hcurl\n", + "\n", + "# The inner product < div v, div E > poses a difficulty for the projection method.\n", + "# In order for us to obtain the correct result, we must modify the transposed gradient\n", + "# as well as the inverse M0 mass matrix\n", + "from psydac.linalg.solvers import inverse\n", + "\n", + "#Gt_0 = P_H1 @ G.T\n", + "\n", + "#M0_0 = P_H1 @ M0 @ P_H1 + P_H1_Gamma\n", + "#M0_inv = inverse(M0_0, 'cg', maxiter=1000, tol=1e-11)\n", + "\n", + "# System Matrix A and A_bc\n", + "#A = # <-\n", + "#A_bc = P_Hcurl @ A @ P_Hcurl + P_Hcurl_Gamma\n", + "\n", + "# rhs vector f and f_bc\n", + "#rho_coeffs = P0(rho_lambdified).coeffs\n", + "#f = # <-\n", + "#f_bc = P_Hcurl @ f" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import time\n", + "\n", + "tol = 1e-10\n", + "maxiter = 1000\n", + "\n", + "# We use a preconditioner for acceleration\n", + "from utils import get_M1_block_kron_solver_2D\n", + "#pc = get_M1_block_kron_solver_2D(V1_h.coeff_space, ncells, degree, periodic=[False, False])\n", + "#M1_0 = P_Hcurl @ M1 @ P_Hcurl + P_Hcurl_Gamma\n", + "#M1_0_inv = inverse(M1_0, 'pcg', pc=pc, tol=1e-11, maxiter=1000)\n", + "#A_bc_inv = inverse(A_bc, 'pcg', pc=M1_0_inv, tol=tol, maxiter=maxiter)\n", + "\n", + "t0 = time.time()\n", + "#E_h = A_bc_inv @ f_bc\n", + "t1 = time.time()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Computing the error norm\n", + "\n", + "As the analytical solution is available, we want to compute the $L^2$ norm of the error.\n", + "In this example, the analytical solution is given by\n", + "\n", + "$$\n", + "\\boldsymbol{E}_{ex}(x, y) = -\\frac{1}{\\pi}\\left(\\begin{matrix} \n", + " \\cos(\\pi x) \\sin(\\pi y) \\\\\n", + " \\cos(\\pi y) \\sin(\\pi x)\n", + " \\end{matrix}\\right)\n", + "$$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from psydac.fem.basic import FemField\n", + "\n", + "from sympde.expr import Norm\n", + "\n", + "E_ex = Tuple( (-1/pi)*cos(pi*x) * sin(pi*y),\n", + " (-1/pi)*cos(pi*y) * sin(pi*x) )\n", + "#E_h_FemField = FemField(V1_h, E_h)\n", + "\n", + "error = Matrix([u1[0] - E_ex[0], u1[1] - E_ex[1]])\n", + "\n", + "# create the formal Norm object\n", + "l2norm = Norm(error, domain, kind='l2')\n", + "\n", + "# discretize the norm\n", + "#l2norm_h = discretize(l2norm, domain_h, V1_h, backend=backend)\n", + "\n", + "# assemble the norm\n", + "#l2_error = l2norm_h.assemble(u1=E_h_FemField)\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( '> CG info :: ',A_bc_inv.get_info() )\n", + "#print( '> L2 error :: {:.2e}'.format( l2_error ) )\n", + "#print( '' )\n", + "#print( '> Solution time :: {:.3g}'.format( t1-t0 ) )" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Visualization\n", + "\n", + "We plot the true solution $\\boldsymbol{E}_{ex}$, the approximate solution $\\boldsymbol{E}_h$ and the error function $|\\boldsymbol{E}_{ex} - \\boldsymbol{E}_h|$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from utils import plot\n", + "\n", + "E_ex_x = lambdify((x, y), E_ex[0])\n", + "E_ex_y = lambdify((x, y), E_ex[1])\n", + "#E_h_x = E_h_FemField[0]\n", + "#E_h_y = E_h_FemField[1]\n", + "#error_x = lambda x, y: abs(E_ex_x(x, y) - E_h_x(x, y))\n", + "#error_y = lambda x, y: abs(E_ex_y(x, y) - E_h_y(x, y))\n", + "\n", + "#plot(gridsize_x = 100, \n", + "# gridsize_y = 100, \n", + "# title = r'Approximation of Solution $\\boldsymbol{E}$, first component', \n", + "# funs = [E_ex_x, E_h_x, error_x], \n", + "# titles = [r'$(\\boldsymbol{E}_{ex})_1(x,y)$', r'$(\\boldsymbol{E}_h)_1(x,y)$', r'$|(\\boldsymbol{E}_{ex}-\\boldsymbol{E}_h)_1(x,y)|$'],\n", + "# surface_plot = True\n", + "#)\n", + "\n", + "#plot(gridsize_x = 100,\n", + "# gridsize_y = 100,\n", + "# title = r'Approximation of Solution $\\boldsymbol{E}$, second component',\n", + "# funs = [E_ex_y, E_h_y, error_y],\n", + "# titles = [r'$(\\boldsymbol{E}_{ex})_2(x,y)$', r'$(\\boldsymbol{E}_h)_2(x,y)$', r'$|(\\boldsymbol{E}_{ex}-\\boldsymbol{E}_h)_2(x,y)|$'],\n", + "# surface_plot = True\n", + "#)" + ] + } + ], + "metadata": {}, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/exercises/harmonic-feec.ipynb b/exercises/harmonic-feec.ipynb new file mode 100644 index 0000000..d7da6a7 --- /dev/null +++ b/exercises/harmonic-feec.ipynb @@ -0,0 +1,270 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "9f28f9af", + "metadata": {}, + "source": [ + "# Harmonic fields with structure-preserving (feec) fem spaces\n", + "\n", + "In this example we .... consider the vector Poisson equation with homogeneous Dirichlet boundary conditions:\n", + "\n", + "$$\n", + "\\begin{align}\n", + " - \\nabla^2 \\mathbf{u} = \\mathbf{f} \\quad \\mbox{in} ~ \\Omega, \\quad \\quad \n", + " \\mathbf{u} = 0 \\quad \\mbox{on} ~ \\partial \\Omega.\n", + "\\end{align}\n", + "$$\n", + "\n", + "## The Variational Formulation\n", + "\n", + "The corresponding variational formulation, using $\\mathbf{H}^1$ formulation, *i.e.* all components are in $H^1$, reads \n", + "\n", + "$$\n", + "\\begin{align}\n", + " \\text{find $\\mathbf{u} \\in V$ such that} \\quad \n", + " a(\\mathbf{u},\\mathbf{v}) = l(\\mathbf{v}) \\quad \\forall \\mathbf{v} \\in V,\n", + "\\end{align}\n", + "$$\n", + "\n", + "where \n", + "\n", + "- $V \\subset \\mathbf{H}_0^1(\\Omega)$, \n", + "- $a(\\mathbf{u},\\mathbf{v}) := \\int_{\\Omega} \\nabla \\mathbf{u} : \\nabla \\mathbf{v} ~ d\\Omega$,\n", + "- $l(\\mathbf{v}) := \\int_{\\Omega} \\mathbf{f} \\cdot \\mathbf{v} ~ d\\Omega$." + ] + }, + { + "cell_type": "markdown", + "id": "5a958607", + "metadata": {}, + "source": [ + "## Formal Model" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d742586c", + "metadata": {}, + "outputs": [], + "source": [ + "from sympde.expr import BilinearForm, LinearForm, integral\n", + "from sympde.expr import find, EssentialBC, Norm, SemiNorm\n", + "from sympde.topology import VectorFunctionSpace, Cube, element_of\n", + "from sympde.calculus import grad, inner, dot\n", + "\n", + "from sympy import pi, sin, Tuple, Matrix\n", + "\n", + "domain = Cube()\n", + "\n", + "V = VectorFunctionSpace('V', domain)\n", + "\n", + "x,y,z = domain.coordinates\n", + "\n", + "u,v = [element_of(V, name=i) for i in ['u', 'v']]\n", + "\n", + "# bilinear form\n", + "a = BilinearForm((u,v), integral(domain , inner(grad(v), grad(u))))\n", + "\n", + "# linear form\n", + "f1 = 3*pi**2*sin(pi*x)*sin(pi*y)*sin(pi*z)\n", + "f2 = 3*pi**2*sin(pi*x)*sin(pi*y)*sin(pi*z)\n", + "f3 = 3*pi**2*sin(pi*x)*sin(pi*y)*sin(pi*z)\n", + "f = Tuple(f1, f2, f3)\n", + "\n", + "l = LinearForm(v, integral(domain, dot(f,v)))\n", + "\n", + "# Dirichlet boundary conditions\n", + "bc = [EssentialBC(u, 0, domain.boundary)]\n", + "\n", + "# Variational problem\n", + "equation = find(u, forall=v, lhs=a(u, v), rhs=l(v), bc=bc)" + ] + }, + { + "cell_type": "markdown", + "id": "4f983ece", + "metadata": {}, + "source": [ + "## Discretization" + ] + }, + { + "cell_type": "markdown", + "id": "51095918", + "metadata": {}, + "source": [ + "We shall need the **discretize** function from **PsyDAC**." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a2a0a2a1", + "metadata": {}, + "outputs": [], + "source": [ + "from psydac.api.discretization import discretize" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "00e54163", + "metadata": {}, + "outputs": [], + "source": [ + "degree = [2,2,2]\n", + "ncells = [8,8,8]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5999c62b", + "metadata": {}, + "outputs": [], + "source": [ + "# Create computational domain from topological domain\n", + "domain_h = discretize(domain, ncells=ncells, comm=None)\n", + "\n", + "# Create discrete spline space\n", + "Vh = discretize(V, domain_h, degree=degree)\n", + "\n", + "# Discretize equation\n", + "equation_h = discretize(equation, domain_h, [Vh, Vh])" + ] + }, + { + "cell_type": "markdown", + "id": "7b29fbcf", + "metadata": {}, + "source": [ + "## Solving the PDE" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "541192ee", + "metadata": {}, + "outputs": [], + "source": [ + "uh = equation_h.solve()" + ] + }, + { + "cell_type": "markdown", + "id": "5174c4b5", + "metadata": {}, + "source": [ + "## Computing the error norm\n", + "\n", + "When the analytical solution is available, you might be interested in computing the $L^2$ norm or $H^1_0$ semi-norm.\n", + "SymPDE allows you to do so, by creating the **Norm** object.\n", + "In this example, the analytical solution is given by\n", + "\n", + "$$\n", + "u_e = \\sin(\\pi x) \\sin(\\pi y) \\sin(\\pi z)\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "3a31c46f", + "metadata": {}, + "source": [ + "### Computing the $L^2$ norm" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5925c6cd", + "metadata": {}, + "outputs": [], + "source": [ + "ue1 = sin(pi*x)*sin(pi*y)*sin(pi*z)\n", + "ue2 = sin(pi*x)*sin(pi*y)*sin(pi*z)\n", + "ue3 = sin(pi*x)*sin(pi*y)*sin(pi*z)\n", + "ue = Tuple(ue1, ue2, ue3)\n", + "\n", + "u = element_of(V, name='u')\n", + "\n", + "error = Matrix([u[0]-ue[0], u[1]-ue[1], u[2]-ue[2]])\n", + "\n", + "# create the formal Norm object\n", + "l2norm = Norm(error, domain, kind='l2')\n", + "\n", + "# discretize the norm\n", + "l2norm_h = discretize(l2norm, domain_h, Vh)\n", + "\n", + "# assemble the norm\n", + "l2_error = l2norm_h.assemble(u=uh)\n", + "\n", + "# print the result\n", + "print(l2_error)" + ] + }, + { + "cell_type": "markdown", + "id": "a6cbfeae", + "metadata": {}, + "source": [ + "### Computing the $H^1$ semi-norm" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e5c1a8b8", + "metadata": {}, + "outputs": [], + "source": [ + "# create the formal Norm object\n", + "h1norm = SemiNorm(error, domain, kind='h1')\n", + "\n", + "# discretize the norm\n", + "h1norm_h = discretize(h1norm, domain_h, Vh)\n", + "\n", + "# assemble the norm\n", + "h1_error = h1norm_h.assemble(u=uh)\n", + "\n", + "# print the result\n", + "print(h1_error)" + ] + }, + { + "cell_type": "markdown", + "id": "3c09131c", + "metadata": {}, + "source": [ + "### Computing the $H^1$ norm" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d829e410", + "metadata": {}, + "outputs": [], + "source": [ + "# create the formal Norm object\n", + "h1norm = Norm(error, domain, kind='h1')\n", + "\n", + "# discretize the norm\n", + "h1norm_h = discretize(h1norm, domain_h, Vh)\n", + "\n", + "# assemble the norm\n", + "h1_error = h1norm_h.assemble(u=uh)\n", + "\n", + "# print the result\n", + "print(h1_error)" + ] + } + ], + "metadata": {}, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/exercises/harmonic-fem.ipynb b/exercises/harmonic-fem.ipynb new file mode 100644 index 0000000..02a87d9 --- /dev/null +++ b/exercises/harmonic-fem.ipynb @@ -0,0 +1,280 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "9f28f9af", + "metadata": {}, + "source": [ + "# Harmonic fields with H1 fem spaces\n", + "\n", + "In this exercise we consider the following vector problem with homogeneous boundary conditions in H(curl):\n", + "\n", + "$$\n", + "\\begin{align}\n", + " \\nabla \\times \\mathbf{u} = 0 \\quad \\mbox{in} ~ \\Omega \\\\\n", + " \\nabla \\cdot \\mathbf{u} = 0 \\quad \\mbox{in} ~ \\Omega \\\\\n", + " \\mathbf{n} \\times \\mathbf{u} = 0 \\quad \\mbox{on} ~ \\partial \\Omega.\n", + "\\end{align}\n", + "$$\n", + "\n", + "## The Variational Formulation\n", + "\n", + "The corresponding variational formulation, using $\\mathbf{H}^1$ formulation, *i.e.* all components are in $H^1$, reads \n", + "\n", + "$$\n", + "\\begin{align}\n", + " \\text{find $\\mathbf{u} \\in V$ such that} \\quad \n", + " a(\\mathbf{u},\\mathbf{v}) = 0 \\quad \\forall \\mathbf{v} \\in V,\n", + "\\end{align}\n", + "$$\n", + "\n", + "where \n", + "\n", + "- $V \\subset (\\mathbf{H}^1(\\Omega))^2$, \n", + "- $a := a_{\\rm curl} + a_{\\rm div} + a_{\\rm bc}$\n", + "- $a_{\\rm curl}(\\mathbf{u},\\mathbf{v}) = \\int_{\\Omega} (\\nabla \\times \\mathbf{u}) \\cdot (\\nabla \\times \\mathbf{v})$,\n", + "- $a_{\\rm div}(\\mathbf{u},\\mathbf{v}) = \\int_{\\Omega} (\\nabla \\cdot \\mathbf{u}) (\\nabla \\cdot \\mathbf{v})$,\n", + "- $a_{\\rm bc}(\\mathbf{u},\\mathbf{v}) = \\kappa \\int_{\\partial\\Omega} (\\mathbf{n} \\times \\mathbf{u}) \\cdot (\\mathbf{n} \\times \\mathbf{v})$.\n", + "\n", + "with $\\kappa$ a large penalization constant.\n", + "\n", + "The exercise is to assemble the matrix of $a$ and compute the first eigenmode using scipy.\n", + "\n", + "The same problem will be addressed in another notebook with structure preserving finite elements, where the solution will only be in H(curl)\n", + "and the divergence-free equation will only be imposed in a weak sense.\n" + ] + }, + { + "cell_type": "markdown", + "id": "5a958607", + "metadata": {}, + "source": [ + "## Formal Model" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d742586c", + "metadata": {}, + "outputs": [], + "source": [ + "from sympde.expr import BilinearForm, LinearForm, integral\n", + "from sympde.expr import find, EssentialBC, Norm, SemiNorm\n", + "from sympde.topology import VectorFunctionSpace, Cube, element_of\n", + "from sympde.calculus import grad, inner, dot\n", + "\n", + "from sympy import pi, sin, Tuple, Matrix\n", + "\n", + "domain = Cube()\n", + "\n", + "V = VectorFunctionSpace('V', domain)\n", + "\n", + "x,y,z = domain.coordinates\n", + "\n", + "u,v = [element_of(V, name=i) for i in ['u', 'v']]\n", + "\n", + "# bilinear form\n", + "a = BilinearForm((u,v), integral(domain , inner(grad(v), grad(u))))\n", + "\n", + "# linear form\n", + "f1 = 3*pi**2*sin(pi*x)*sin(pi*y)*sin(pi*z)\n", + "f2 = 3*pi**2*sin(pi*x)*sin(pi*y)*sin(pi*z)\n", + "f3 = 3*pi**2*sin(pi*x)*sin(pi*y)*sin(pi*z)\n", + "f = Tuple(f1, f2, f3)\n", + "\n", + "l = LinearForm(v, integral(domain, dot(f,v)))\n", + "\n", + "# Dirichlet boundary conditions\n", + "bc = [EssentialBC(u, 0, domain.boundary)]\n", + "\n", + "# Variational problem\n", + "equation = find(u, forall=v, lhs=a(u, v), rhs=l(v), bc=bc)" + ] + }, + { + "cell_type": "markdown", + "id": "4f983ece", + "metadata": {}, + "source": [ + "## Discretization" + ] + }, + { + "cell_type": "markdown", + "id": "51095918", + "metadata": {}, + "source": [ + "We shall need the **discretize** function from **PsyDAC**." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a2a0a2a1", + "metadata": {}, + "outputs": [], + "source": [ + "from psydac.api.discretization import discretize" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "00e54163", + "metadata": {}, + "outputs": [], + "source": [ + "degree = [2,2,2]\n", + "ncells = [8,8,8]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5999c62b", + "metadata": {}, + "outputs": [], + "source": [ + "# Create computational domain from topological domain\n", + "domain_h = discretize(domain, ncells=ncells, comm=None)\n", + "\n", + "# Create discrete spline space\n", + "Vh = discretize(V, domain_h, degree=degree)\n", + "\n", + "# Discretize equation\n", + "equation_h = discretize(equation, domain_h, [Vh, Vh])" + ] + }, + { + "cell_type": "markdown", + "id": "7b29fbcf", + "metadata": {}, + "source": [ + "## Solving the PDE" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "541192ee", + "metadata": {}, + "outputs": [], + "source": [ + "uh = equation_h.solve()" + ] + }, + { + "cell_type": "markdown", + "id": "5174c4b5", + "metadata": {}, + "source": [ + "## Computing the error norm\n", + "\n", + "When the analytical solution is available, you might be interested in computing the $L^2$ norm or $H^1_0$ semi-norm.\n", + "SymPDE allows you to do so, by creating the **Norm** object.\n", + "In this example, the analytical solution is given by\n", + "\n", + "$$\n", + "u_e = \\sin(\\pi x) \\sin(\\pi y) \\sin(\\pi z)\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "3a31c46f", + "metadata": {}, + "source": [ + "### Computing the $L^2$ norm" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5925c6cd", + "metadata": {}, + "outputs": [], + "source": [ + "ue1 = sin(pi*x)*sin(pi*y)*sin(pi*z)\n", + "ue2 = sin(pi*x)*sin(pi*y)*sin(pi*z)\n", + "ue3 = sin(pi*x)*sin(pi*y)*sin(pi*z)\n", + "ue = Tuple(ue1, ue2, ue3)\n", + "\n", + "u = element_of(V, name='u')\n", + "\n", + "error = Matrix([u[0]-ue[0], u[1]-ue[1], u[2]-ue[2]])\n", + "\n", + "# create the formal Norm object\n", + "l2norm = Norm(error, domain, kind='l2')\n", + "\n", + "# discretize the norm\n", + "l2norm_h = discretize(l2norm, domain_h, Vh)\n", + "\n", + "# assemble the norm\n", + "l2_error = l2norm_h.assemble(u=uh)\n", + "\n", + "# print the result\n", + "print(l2_error)" + ] + }, + { + "cell_type": "markdown", + "id": "a6cbfeae", + "metadata": {}, + "source": [ + "### Computing the $H^1$ semi-norm" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e5c1a8b8", + "metadata": {}, + "outputs": [], + "source": [ + "# create the formal Norm object\n", + "h1norm = SemiNorm(error, domain, kind='h1')\n", + "\n", + "# discretize the norm\n", + "h1norm_h = discretize(h1norm, domain_h, Vh)\n", + "\n", + "# assemble the norm\n", + "h1_error = h1norm_h.assemble(u=uh)\n", + "\n", + "# print the result\n", + "print(h1_error)" + ] + }, + { + "cell_type": "markdown", + "id": "3c09131c", + "metadata": {}, + "source": [ + "### Computing the $H^1$ norm" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d829e410", + "metadata": {}, + "outputs": [], + "source": [ + "# create the formal Norm object\n", + "h1norm = Norm(error, domain, kind='h1')\n", + "\n", + "# discretize the norm\n", + "h1norm_h = discretize(h1norm, domain_h, Vh)\n", + "\n", + "# assemble the norm\n", + "h1_error = h1norm_h.assemble(u=uh)\n", + "\n", + "# print the result\n", + "print(h1_error)" + ] + } + ], + "metadata": {}, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/exercises/hodge_laplace_annulus.ipynb b/exercises/hodge_laplace_annulus.ipynb new file mode 100644 index 0000000..2009457 --- /dev/null +++ b/exercises/hodge_laplace_annulus.ipynb @@ -0,0 +1,539 @@ +{ + "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": {}, + "source": [ + "## Obtain the Annulus Domain" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "\n", + "from sympde.topology import Square, PolarMapping\n", + "from sympde.utilities.utils import plot_domain\n", + "\n", + "# rmin and rmax define the inner and outer radius of the annulus\n", + "rmin, rmax = 0.3, 1.\n", + "\n", + "logical_domain = Square('S', bounds1=(0., 1.), bounds2=(0., 2*np.pi))\n", + "F = PolarMapping('A', dim=2, c1=0., c2=0., rmin=rmin, rmax=rmax)\n", + "domain = F(logical_domain)\n", + "\n", + "# periodicity of the domain\n", + "periodic = [False, True] \n", + "\n", + "plot_domain(domain, draw=True, isolines=True)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Discrete Model using de Rham objects\n", + "\n", + "### Derham and Mass Matrices First" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "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", + "\n", + "V0 = derham.V0\n", + "V1 = derham.V1\n", + "V2 = derham.V2\n", + "\n", + "u0, v0 = elements_of(V0, names='u0, v0')\n", + "u1, v1 = elements_of(V1, names='u1, v1')\n", + "u2, v2 = elements_of(V2, names='u2, v2')\n", + "\n", + "# bilinear forms corresponding to V0, V1 and V2 mass matrices\n", + "m0 = BilinearForm((u0, v0), integral(domain, u0*v0))\n", + "m1 = BilinearForm((u1, v1), integral(domain, dot(u1, v1)))\n", + "m2 = BilinearForm((u2, v2), integral(domain, u2*v2))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from psydac.api.discretization import discretize\n", + "from psydac.api.settings import PSYDAC_BACKEND_GPYCCEL\n", + "\n", + "backend = PSYDAC_BACKEND_GPYCCEL" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ncells = [32, 32] # Bspline cells\n", + "degree = [4, 4] # Bspline degree" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# discretize domain and derham\n", + "domain_h = discretize(domain, ncells=ncells, periodic=periodic)\n", + "derham_h = discretize(derham, domain_h, degree=degree)\n", + "\n", + "# define FEM spaces V0_h, V1_h and V2_h\n", + "V0h = derham_h.V0\n", + "V1h = derham_h.V1\n", + "V2h = derham_h.V2\n", + "\n", + "# Commuting projection operators\n", + "P0, P1, P2 = derham_h.projectors()\n", + "\n", + "# Exterior derivative operators (grad and curl)\n", + "G, C = derham_h.derivatives_as_matrices\n", + "Gt = G.T\n", + "Ct = C.T\n", + "\n", + "# Mass matrices\n", + "m0_h = discretize(m0, domain_h, (V0h, V0h), backend=backend)\n", + "m1_h = discretize(m1, domain_h, (V1h, V1h), backend=backend)\n", + "m2_h = discretize(m2, domain_h, (V2h, V2h), backend=backend)\n", + "\n", + "M0 = m0_h.assemble()\n", + "M1 = m1_h.assemble()\n", + "M2 = m2_h.assemble()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Boundary Conditions\n", + "\n", + "We choose to apply the projection method. For that matter, we use the projection matrices\n", + "$$\n", + "\\begin{align*}\n", + " &\\mathbb{P}_{H_0^1}, \\text{ projecting into the coefficients space of functions in }H_0^1 \\\\\n", + " &\\mathbb{P}_{H_0(curl)}, \\text{ projecting into the coefficients space of functions in }H_0(curl)\n", + "\\end{align*}\n", + "$$\n", + "as well as the projection matrices\n", + "$$\n", + "\\begin{align*}\n", + " &\\mathbb{P}_{H_0^1}^{\\Gamma} = \\mathbb{I}_0 - \\mathbb{P}_{H_0^1} \\\\\n", + " &\\mathbb{P}_{H_0(curl)}^{\\Gamma} = \\mathbb{I}_1 - \\mathbb{P}_{H_0(curl)}\n", + "\\end{align*}\n", + "$$\n", + "The latter two matrices are used for stabilization purposes, i.e., in order to keep the system matrix non-singular." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from psydac.linalg.basic import IdentityOperator\n", + "\n", + "from utils import H1BoundaryProjector2D, HcurlBoundaryProjector2D\n", + "\n", + "P_H1 = H1BoundaryProjector2D(V0, V0h.coeff_space, periodic=periodic)\n", + "P_Hcurl = HcurlBoundaryProjector2D(V1, V1h.coeff_space, periodic=periodic)\n", + "\n", + "I0 = IdentityOperator(V0h.coeff_space)\n", + "I1 = IdentityOperator(V1h.coeff_space)\n", + "\n", + "P_H1_Gamma = I0 - P_H1\n", + "P_Hcurl_Gamma = I1 - P_Hcurl" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 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", + "\\begin{align*}\n", + " \\mathbb{M}_{I,I} x_{I} = b_{I}\n", + "\\end{align*}\n", + "$$\n", + "but we only have access to the larger matrices and vectors\n", + "$$\n", + "\\begin{align*}\n", + " &\\mathbb{M} = \\left(\\begin{matrix} \\mathbb{M}_{I,I} && \\mathbb{M}_{I,\\Gamma} \\\\\n", + " \\mathbb{M}_{\\Gamma,I} && \\mathbb{M}_{\\Gamma,\\Gamma}\n", + " \\end{matrix}\\right), \\\\\n", + " &b = \\left(\\begin{matrix} b_{I} \\\\\n", + " b_{\\Gamma}\n", + " \\end{matrix}\\right).\n", + "\\end{align*}\n", + "$$\n", + "One way to do this is to solve the following system\n", + "$$\n", + "\\begin{align*}\n", + " \\left(\\begin{matrix} \\mathbb{M}_{I,I} && \\mathbb{0} \\\\\n", + " \\mathbb{0} && \\mathbb{I}_{\\Gamma,\\Gamma}\n", + " \\end{matrix}\\right)\n", + " \\left(\\begin{matrix} x_{I} \\\\ x_{\\Gamma} \\end{matrix}\\right) = \n", + " \\left(\\begin{matrix} b_{I} \\\\ 0 \\end{matrix}\\right)\n", + "\\end{align*}\n", + "$$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "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-12)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Build the discrete $\\mathcal{L}_h^1$ operator in sparse format\n", + "\n", + "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", + "\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", + "\\mathcal{L}_h^1 = \\widetilde{\\boldsymbol{curl}}_h\\ curl - \\boldsymbol{grad}\\ \\widetilde{div}_h\n", + "$$\n", + "in order to later enforce the orthogonality constraint." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from scipy.sparse import csc_matrix, bmat\n", + "from scipy.sparse.linalg import inv\n", + "\n", + "# Sparse representations of projection operators\n", + "P_H1_sp = csc_matrix(P_H1.toarray())\n", + "P_Hcurl_sp = csc_matrix(P_Hcurl.toarray())\n", + "\n", + "I0_sp = csc_matrix(I0.toarray())\n", + "I1_sp = csc_matrix(I1.toarray())\n", + "\n", + "P_H1_Gamma_sp = I0_sp - P_H1_sp\n", + "P_Hcurl_Gamma_sp = I1_sp - P_Hcurl_sp\n", + "\n", + "# Sparse representations of mass matrices and differentiation operators\n", + "M0_sp = M0.tosparse()\n", + "M1_sp = M1.tosparse()\n", + "M2_sp = M2.tosparse()\n", + "\n", + "G_sp = G.tosparse()\n", + "C_sp = C.tosparse()\n", + "\n", + "Gt_sp = G_sp.transpose()\n", + "Ct_sp = C_sp.transpose()\n", + "\n", + "# It can be shown that for the projection method to work \n", + "# we need a modified version of the gradient operator\n", + "G_sp_0 = G.tosparse() @ P_H1_sp\n", + "Gt_sp_0 = G_sp_0.transpose()\n", + "\n", + "# ... and of course the correct inverse M0 mass matrix in its sparse representation\n", + "M0_sp_0 = P_H1_sp @ M0.tosparse() @ P_H1_sp + P_H1_Gamma_sp\n", + "inv_M0_sp = inv(M0_sp_0.tocsc())\n", + "inv_M0_sp.eliminate_zeros()\n", + "\n", + "# For this complex operator, the projection method consists of using the modified gradient\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" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Compute first eigenvector of $\\mathcal{L}_h^1$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "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", + "\n", + "# A basis of the vector space of harmonic forms (hf)\n", + "#hf = eigenvectors[:,0]\n", + "# Obtain the eigenvector as a \"row csc_matrix (1 x n)\"\n", + "#hft_sp = csc_matrix(hf)\n", + "# ... and as a \"column csc_matrix (n x 1)\"\n", + "#hf_sp = hft_sp.transpose()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Build system matrix in sparse format\n", + "... and apply projection method to enforce boundary conditions" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# System matrix without taking BCs into account\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", + "\n", + "P = bmat([[P_H1_sp, None, None],\n", + " [None, P_Hcurl_sp, None],\n", + " [None, None, ONE]])\n", + "\n", + "P_Gamma = bmat([[P_H1_Gamma_sp, None, None],\n", + " [None, P_Hcurl_Gamma_sp, None],\n", + " [None, None, ZERO]])\n", + "\n", + "# System matrix taking BCs into account\n", + "#A_sp_bc = P @ A_sp @ P + P_Gamma" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "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", + "#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 @ 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", + "\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", + "t1 = time.time()\n", + "# --------------------------------------------------------------" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Error Analysis" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# 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", + "#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", + "from psydac.fem.plotting_utilities import plot_field_2d\n", + "\n", + "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')" + ] + } + ], + "metadata": {}, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/exercises/poisson_2d_I.ipynb b/exercises/poisson_2d_I.ipynb new file mode 100644 index 0000000..ca64bab --- /dev/null +++ b/exercises/poisson_2d_I.ipynb @@ -0,0 +1,322 @@ +{ + "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 Poisson Solver on the Unit Square\n", + "## ... using Psydac's bilinear form interface\n", + "\n", + "In this exercise we write a solver for the 2D Poisson problem on the unit square using Psydac's bilinear form interface.\n", + "\n", + "\\begin{align*}\n", + " -\\Delta\\phi &= f \\quad \\text{ in }\\Omega=]0,1[^2, \\\\\n", + " \\phi &= 0 \\quad \\text{ on }\\partial\\Omega,\n", + "\\end{align*}\n", + "for the specific choice\n", + "\\begin{align*}\n", + " f(x, y) = [ (-x^4 + x^3 + (-y^2 + \\pi^2)x^2 + (y^2 - 4y - \\pi^2)x + 2y - 2) \\sin(\\pi y) + (2\\pi x^2(1-x)) \\cos(\\pi y)] e^{xy}\n", + "\\end{align*}\n", + "\n", + "## The Variational Formulation\n", + "\n", + "The corresponding variational formulation reads\n", + "\n", + "\\begin{align*}\n", + " \\text{Find }\\phi\\in V\\coloneqq H_0^1(\\Omega)\\text{ such that } \\\\\n", + " a(\\phi, \\psi) = l(\\psi)\\quad\\forall\\ \\psi\\in V\n", + "\\end{align*}\n", + "\n", + "where \n", + "- $a(\\phi,\\psi) \\coloneqq \\int_{\\Omega} \\nabla\\phi\\cdot\\nabla\\psi ~ d\\Omega$,\n", + "- $l(\\psi) := \\int_{\\Omega} f\\psi ~ d\\Omega$." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Formal Model" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from sympde.calculus import dot, grad\n", + "from sympde.expr import BilinearForm, LinearForm, integral\n", + "from sympde.topology import elements_of, Square, ScalarFunctionSpace\n", + "\n", + "domain = Square('S', bounds1=(0, 1), bounds2=(0, 1))\n", + "\n", + "V = ScalarFunctionSpace('V', domain, kind='h1')\n", + "\n", + "phi, psi = elements_of(V, names='phi, psi')\n", + "\n", + "# bilinear form\n", + "a = BilinearForm((phi, psi), integral(domain, phi*psi))\n", + "\n", + "# linear form\n", + "from sympy import pi, sin, cos, exp\n", + "x,y = domain.coordinates\n", + "f = ( (-x**4 + x**3 + (-y**2 + pi**2)*x**2 + (y**2 - 4*y - pi**2)*x + 2*y - 2) * sin(pi*y) + (2*pi*x**2*(1-x)) * cos(pi*y) ) * exp(x*y)\n", + "\n", + "l = LinearForm(psi, integral(domain, f*psi))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Discretization\n", + "\n", + "We will use Psydac to discretize the problem." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from psydac.api.discretization import discretize\n", + "from psydac.api.settings import PSYDAC_BACKEND_GPYCCEL\n", + "\n", + "backend = PSYDAC_BACKEND_GPYCCEL" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ncells = [16, 16] # Bspline cells\n", + "degree = [3, 3] # Bspline degree" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "domain_h = discretize(domain, ncells=ncells, periodic=[False, False])\n", + "V_h = discretize(V, domain_h, degree=degree)\n", + "\n", + "a_h = discretize(a, domain_h, (V_h, V_h), backend=backend)\n", + "l_h = discretize(l, domain_h, V_h, backend=backend)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Boundary Conditions\n", + "\n", + "We choose to apply a projection method. For that matter, we construct the projection matrix $\\mathbb{P}_0$ and its counterpart $\\mathbb{P}_{\\Gamma} = \\mathbb{I} - \\mathbb{P}_0$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from utils import H1BoundaryProjector2D\n", + "from psydac.linalg.basic import IdentityOperator\n", + "\n", + "P_0 = H1BoundaryProjector2D(V, V_h.coeff_space)\n", + "\n", + "I_0 = IdentityOperator(V_h.coeff_space)\n", + "P_Gamma = I_0 - P_0" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "A = a_h.assemble()\n", + "f = l_h.assemble()\n", + "\n", + "A_bc = P_0 @ A @ P_0 + P_Gamma\n", + "f_bc = P_0 @ f" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Solving the PDE" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import time\n", + "\n", + "from psydac.linalg.solvers import inverse\n", + "\n", + "tol = 1e-12\n", + "maxiter = 1000\n", + "\n", + "A_bc_inv = inverse(A_bc, 'cg', tol=tol, maxiter=maxiter)\n", + "\n", + "t0 = time.time()\n", + "phi_h = A_bc_inv @ f_bc\n", + "t1 = time.time()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Computing the error norm\n", + "\n", + "When the analytical solution is available, you might be interested in computing the $L^2$ norm or $H^1_0$ semi-norm.\n", + "SymPDE allows you to do so, by creating the **Norm** object.\n", + "In this example, the analytical solution is given by\n", + "\n", + "$$\n", + "phi_{ex}(x, y) = x(x-1)\\sin(\\pi y)e^{xy}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Computing the $L^2$ norm" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from psydac.fem.basic import FemField\n", + "\n", + "from sympde.expr import Norm, SemiNorm\n", + "\n", + "phi_ex = x*(x-1)*sin(pi*y)*exp(x*y)\n", + "phi_h_FemField = FemField(V_h, phi_h)\n", + "\n", + "error = phi_ex - phi\n", + "\n", + "# create the formal Norm object\n", + "l2norm = Norm(error, domain, kind='l2')\n", + "\n", + "# discretize the norm\n", + "l2norm_h = discretize(l2norm, domain_h, V_h, backend=backend)\n", + "\n", + "# assemble the norm\n", + "l2_error = l2norm_h.assemble(phi=phi_h_FemField)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Computing the $H^1$ semi-norm" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# create the formal Norm object\n", + "h1norm = SemiNorm(error, domain, kind='h1')\n", + "\n", + "# discretize the norm\n", + "h1norm_h = discretize(h1norm, domain_h, V_h)\n", + "\n", + "# assemble the norm\n", + "h1semi_error = h1norm_h.assemble(phi=phi_h_FemField)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Computing the $H^1$ norm" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# create the formal Norm object\n", + "h1norm = Norm(error, domain, kind='h1')\n", + "\n", + "# discretize the norm\n", + "h1norm_h = discretize(h1norm, domain_h, V_h)\n", + "\n", + "# assemble the norm\n", + "h1_error = h1norm_h.assemble(phi=phi_h_FemField)\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( '> CG info :: ',A_bc_inv.get_info() )\n", + "print( '> L2 error :: {:.2e}'.format( l2_error ) )\n", + "print( '> H1-Semi error :: {:.2e}'.format( h1semi_error ) )\n", + "print( '> H1 error :: {:.2e}'.format( h1_error ) )\n", + "print( '' )\n", + "print( '> Solution time :: {:.3g}'.format( t1-t0 ) )" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Visualization\n", + "\n", + "We plot the true solution $\\phi_{ex}$, the approximate solution $\\phi_h$ and the error function $|\\phi_{ex} - \\phi_h|$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from sympy import lambdify\n", + "\n", + "from utils import plot\n", + "\n", + "phi_ex_fun = lambdify(domain.coordinates, phi_ex)\n", + "error = lambda x, y: abs(phi_ex_fun(x, y) - phi_h_FemField(x, y))\n", + "\n", + "plot(gridsize_x = 100, \n", + " gridsize_y = 100, \n", + " title = r'Approximation of Solution $\\phi$', \n", + " funs = [phi_ex_fun, phi_h_FemField, error], \n", + " titles = [r'$\\phi_{ex}(x,y)$', r'$\\phi_h(x,y)$', r'$|(\\phi_{ex}-\\phi_h)(x,y)|$'],\n", + " surface_plot = True\n", + ")" + ] + } + ], + "metadata": {}, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/exercises/poisson_2d_II.ipynb b/exercises/poisson_2d_II.ipynb new file mode 100644 index 0000000..e22c2d8 --- /dev/null +++ b/exercises/poisson_2d_II.ipynb @@ -0,0 +1,329 @@ +{ + "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 Poisson Solver on the Unit Square\n", + "## ... using Psydac's de Rham interface\n", + "\n", + "In this exercise we write a solver for the 2D Poisson problem on the unit square using Psydac's bilinear form interface.\n", + "\n", + "\\begin{align*}\n", + " -\\Delta\\phi &= f \\quad \\text{ in }\\Omega=]0,1[^2, \\\\\n", + " \\phi &= 0 \\quad \\text{ on }\\partial\\Omega,\n", + "\\end{align*}\n", + "for the specific choice\n", + "\\begin{align*}\n", + " f(x, y) = [ (-x^4 + x^3 + (-y^2 + \\pi^2)x^2 + (y^2 - 4y - \\pi^2)x + 2y - 2) \\sin(\\pi y) + (2\\pi x^2(1-x)) \\cos(\\pi y)] e^{xy}\n", + "\\end{align*}\n", + "\n", + "## The Variational Formulation\n", + "\n", + "The corresponding variational formulation reads\n", + "\n", + "\\begin{align*}\n", + " \\text{Find }\\phi\\in V\\coloneqq H_0^1(\\Omega)\\text{ such that } \\\\\n", + " a(\\phi, \\psi) = l(\\psi)\\quad\\forall\\ \\psi\\in V\n", + "\\end{align*}\n", + "\n", + "where \n", + "- $a(\\phi,\\psi) \\coloneqq \\int_{\\Omega} \\nabla\\phi\\cdot\\nabla\\psi ~ d\\Omega$,\n", + "- $l(\\psi) := \\int_{\\Omega} f\\psi ~ d\\Omega$." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Discrete Model using de Rham objects" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from sympde.calculus import dot, grad\n", + "from sympde.expr import BilinearForm, LinearForm, integral\n", + "from sympde.topology import elements_of, Square, ScalarFunctionSpace, Derham\n", + "\n", + "domain = Square('S', bounds1=(0, 1), bounds2=(0, 1))\n", + "derham = Derham(domain, sequence=['h1', 'hcurl', 'l2'])\n", + "\n", + "V0 = derham.V0\n", + "V1 = derham.V1\n", + "\n", + "u0, v0 = elements_of(V0, names='u0, v0')\n", + "u1, v1 = elements_of(V1, names='u1, v1')\n", + "\n", + "# bilinear form corresponding to V1 mass matrix\n", + "m1 = BilinearForm((u1, v1), integral(domain, dot(u1, v1)))\n", + "\n", + "# linear form\n", + "from sympy import pi, sin, cos, exp\n", + "x,y = domain.coordinates\n", + "f = ( (-x**4 + x**3 + (-y**2 + pi**2)*x**2 + (y**2 - 4*y - pi**2)*x + 2*y - 2) * sin(pi*y) + (2*pi*x**2*(1-x)) * cos(pi*y) ) * exp(x*y)\n", + "\n", + "l = LinearForm(v0, integral(domain, f*v0))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from psydac.api.discretization import discretize\n", + "from psydac.api.settings import PSYDAC_BACKEND_GPYCCEL\n", + "\n", + "backend = PSYDAC_BACKEND_GPYCCEL" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ncells = [16, 16] # Bspline cells\n", + "degree = [3, 3] # Bspline degree" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "domain_h = discretize(domain, ncells=ncells, periodic=[False, False])\n", + "derham_h = discretize(derham, domain_h, degree=degree)\n", + "\n", + "V0_h = derham_h.V0\n", + "V1_h = derham_h.V1\n", + "\n", + "# Exterior Derivative operators (grad and curl)\n", + "G, C = derham_h.derivatives_as_matrices\n", + "\n", + "# Mass Matrix\n", + "m1_h = discretize(m1, domain_h, (V1_h, V1_h), backend=backend)\n", + "\n", + "M1 = m1_h.assemble()\n", + "\n", + "# System Matrix A\n", + "from psydac.linalg.basic import IdentityOperator\n", + "A = G.T @ M1 @ G\n", + "\n", + "# rhs vector f\n", + "l_h = discretize(l, domain_h, V0_h, backend=backend)\n", + "f = l_h.assemble()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Boundary Conditions\n", + "\n", + "We choose to apply a projection method. For that matter, we construct the projection matrix $\\mathbb{P}_0$ and its counterpart $\\mathbb{P}_{\\Gamma} = \\mathbb{I} - \\mathbb{P}_0$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from utils import H1BoundaryProjector2D\n", + "from psydac.linalg.basic import IdentityOperator\n", + "\n", + "P_0 = H1BoundaryProjector2D(V0, V0_h.coeff_space)\n", + "\n", + "I_0 = IdentityOperator(V0_h.coeff_space)\n", + "P_Gamma = I_0 - P_0" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "A_bc = P_0 @ A @ P_0 + P_Gamma\n", + "f_bc = P_0 @ f" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Solving the PDE" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import time\n", + "\n", + "from psydac.linalg.solvers import inverse\n", + "\n", + "tol = 1e-12\n", + "maxiter = 1000\n", + "\n", + "A_bc_inv = inverse(A_bc, 'cg', tol=tol, maxiter=maxiter)\n", + "\n", + "t0 = time.time()\n", + "phi_h = A_bc_inv @ f_bc\n", + "t1 = time.time()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Computing the error norm\n", + "\n", + "When the analytical solution is available, you might be interested in computing the $L^2$ norm or $H^1_0$ semi-norm.\n", + "SymPDE allows you to do so, by creating the **Norm** object.\n", + "In this example, the analytical solution is given by\n", + "\n", + "$$\n", + "phi_{ex}(x, y) = x(x-1)\\sin(\\pi y)e^{xy}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Computing the $L^2$ norm" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from psydac.fem.basic import FemField\n", + "\n", + "from sympde.expr import Norm, SemiNorm\n", + "\n", + "phi_ex = x*(x-1)*sin(pi*y)*exp(x*y)\n", + "phi_h_FemField = FemField(V0_h, phi_h)\n", + "\n", + "error = phi_ex - u0\n", + "\n", + "# create the formal Norm object\n", + "l2norm = Norm(error, domain, kind='l2')\n", + "\n", + "# discretize the norm\n", + "l2norm_h = discretize(l2norm, domain_h, V0_h, backend=backend)\n", + "\n", + "# assemble the norm\n", + "l2_error = l2norm_h.assemble(u0=phi_h_FemField)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Computing the $H^1$ semi-norm" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# create the formal Norm object\n", + "h1norm = SemiNorm(error, domain, kind='h1')\n", + "\n", + "# discretize the norm\n", + "h1norm_h = discretize(h1norm, domain_h, V0_h)\n", + "\n", + "# assemble the norm\n", + "h1semi_error = h1norm_h.assemble(u0=phi_h_FemField)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Computing the $H^1$ norm" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# create the formal Norm object\n", + "h1norm = Norm(error, domain, kind='h1')\n", + "\n", + "# discretize the norm\n", + "h1norm_h = discretize(h1norm, domain_h, V0_h)\n", + "\n", + "# assemble the norm\n", + "h1_error = h1norm_h.assemble(u0=phi_h_FemField)\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( '> CG info :: ',A_bc_inv.get_info() )\n", + "print( '> L2 error :: {:.2e}'.format( l2_error ) )\n", + "print( '> H1-Semi error :: {:.2e}'.format( h1semi_error ) )\n", + "print( '> H1 error :: {:.2e}'.format( h1_error ) )\n", + "print( '' )\n", + "print( '> Solution time :: {:.3g}'.format( t1-t0 ) )" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Visualization\n", + "\n", + "We plot the true solution $\\phi_{ex}$, the approximate solution $\\phi_h$ and the error function $|\\phi_{ex} - \\phi_h|$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from sympy import lambdify\n", + "\n", + "from utils import plot\n", + "\n", + "phi_ex_fun = lambdify(domain.coordinates, phi_ex)\n", + "error = lambda x, y: abs(phi_ex_fun(x, y) - phi_h_FemField(x, y))\n", + "\n", + "plot(gridsize_x = 100, \n", + " gridsize_y = 100, \n", + " title = r'Approximation of Solution $\\phi$', \n", + " funs = [phi_ex_fun, phi_h_FemField, error], \n", + " titles = [r'$\\phi_{ex}(x,y)$', r'$\\phi_h(x,y)$', r'$|(\\phi_{ex}-\\phi_h)(x,y)|$'],\n", + " surface_plot = True\n", + ")" + ] + } + ], + "metadata": {}, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/exercises/poisson_2d_I_naturalBC.ipynb b/exercises/poisson_2d_I_naturalBC.ipynb new file mode 100644 index 0000000..58e6866 --- /dev/null +++ b/exercises/poisson_2d_I_naturalBC.ipynb @@ -0,0 +1,297 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# 2D Poisson Solver on the Unit Square\n", + "## ... using Psydac's bilinear form interface\n", + "\n", + "In this exercise we write a solver for the 2D Poisson problem with natural BCs on the unit square using Psydac's bilinear form interface.\n", + "\n", + "\\begin{align*}\n", + " -\\Delta\\phi &= f \\quad \\text{ in }\\Omega=]0,1[^2, \\\\\n", + " \\frac{\\partial\\phi}{\\partial \\boldsymbol{n}} &= 0 \\quad \\text{ on }\\partial\\Omega,\n", + "\\end{align*}\n", + "for the specific choice\n", + "\\begin{align*}\n", + " f(x, y) = 8\\pi^2\\cos(2\\pi x)\\cos(2\\pi y)\n", + "\\end{align*}\n", + "\n", + "## The Variational Formulation\n", + "\n", + "The corresponding variational formulation reads\n", + "\n", + "\\begin{align*}\n", + " \\text{Find }\\phi\\in V\\coloneqq H^1(\\Omega)\\text{ such that } \\\\\n", + " a(\\phi, \\psi) = l(\\psi)\\quad\\forall\\ \\psi\\in V\n", + "\\end{align*}\n", + "\n", + "where \n", + "- $a(\\phi,\\psi) \\coloneqq \\int_{\\Omega} \\nabla\\phi\\cdot\\nabla\\psi ~ d\\Omega$,\n", + "- $l(\\psi) := \\int_{\\Omega} f\\psi ~ d\\Omega$." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Formal Model" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from sympde.calculus import dot, grad\n", + "from sympde.expr import BilinearForm, LinearForm, integral\n", + "from sympde.topology import elements_of, Square, ScalarFunctionSpace\n", + "\n", + "domain = Square('S', bounds1=(0, 1), bounds2=(0, 1))\n", + "\n", + "V = ScalarFunctionSpace('V', domain, kind='h1')\n", + "\n", + "phi, psi = elements_of(V, names='phi, psi')\n", + "\n", + "# bilinear form\n", + "a = BilinearForm((phi, psi), integral(domain, dot(grad(phi), grad(psi))))\n", + "\n", + "# linear form\n", + "from sympy import pi, sin, cos, exp\n", + "x,y = domain.coordinates\n", + "f = 8 * pi**2 * cos(2*pi*x) * cos(2*pi*y)\n", + "\n", + "l = LinearForm(psi, integral(domain, f*psi))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Discretization\n", + "\n", + "We will use Psydac to discretize the problem." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from psydac.api.discretization import discretize\n", + "from psydac.api.settings import PSYDAC_BACKEND_GPYCCEL\n", + "\n", + "backend = PSYDAC_BACKEND_GPYCCEL" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ncells = [16, 16] # Bspline cells\n", + "degree = [3, 3] # Bspline degree" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "domain_h = discretize(domain, ncells=ncells, periodic=[False, False])\n", + "V_h = discretize(V, domain_h, degree=degree)\n", + "\n", + "a_h = discretize(a, domain_h, (V_h, V_h), backend=backend)\n", + "l_h = discretize(l, domain_h, V_h, backend=backend)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Boundary Conditions\n", + "\n", + "We must not take care of boundary conditions. They are included **naturally** in the variational formulation." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "A = a_h.assemble()\n", + "f = l_h.assemble()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Solving the PDE" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import time\n", + "\n", + "from psydac.linalg.solvers import inverse\n", + "\n", + "tol = 1e-12\n", + "maxiter = 1000\n", + "\n", + "A_bc_inv = inverse(A, 'cg', tol=tol, maxiter=maxiter)\n", + "\n", + "t0 = time.time()\n", + "phi_h = A_bc_inv @ f\n", + "t1 = time.time()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Computing the error norm\n", + "\n", + "When the analytical solution is available, you might be interested in computing the $L^2$ norm or $H^1_0$ semi-norm.\n", + "SymPDE allows you to do so, by creating the **Norm** object.\n", + "In this example, the analytical solution is given by\n", + "\n", + "$$\n", + "phi_{ex}(x, y) = \\cos(2\\pi x)\\cos(2\\pi y)\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Computing the $L^2$ norm" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from psydac.fem.basic import FemField\n", + "\n", + "from sympde.expr import Norm, SemiNorm\n", + "\n", + "phi_ex = cos(2*pi*x) * cos(2*pi*y)\n", + "phi_h_FemField = FemField(V_h, phi_h)\n", + "\n", + "error = phi_ex - phi\n", + "\n", + "# create the formal Norm object\n", + "l2norm = Norm(error, domain, kind='l2')\n", + "\n", + "# discretize the norm\n", + "l2norm_h = discretize(l2norm, domain_h, V_h, backend=backend)\n", + "\n", + "# assemble the norm\n", + "l2_error = l2norm_h.assemble(phi=phi_h_FemField)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Computing the $H^1$ semi-norm" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# create the formal Norm object\n", + "h1norm = SemiNorm(error, domain, kind='h1')\n", + "\n", + "# discretize the norm\n", + "h1norm_h = discretize(h1norm, domain_h, V_h)\n", + "\n", + "# assemble the norm\n", + "h1semi_error = h1norm_h.assemble(phi=phi_h_FemField)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Computing the $H^1$ norm" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# create the formal Norm object\n", + "h1norm = Norm(error, domain, kind='h1')\n", + "\n", + "# discretize the norm\n", + "h1norm_h = discretize(h1norm, domain_h, V_h)\n", + "\n", + "# assemble the norm\n", + "h1_error = h1norm_h.assemble(phi=phi_h_FemField)\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( '> CG info :: ',A_bc_inv.get_info() )\n", + "print( '> L2 error :: {:.2e}'.format( l2_error ) )\n", + "print( '> H1-Semi error :: {:.2e}'.format( h1semi_error ) )\n", + "print( '> H1 error :: {:.2e}'.format( h1_error ) )\n", + "print( '' )\n", + "print( '> Solution time :: {:.3g}'.format( t1-t0 ) )" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Visualization\n", + "\n", + "We plot the true solution $\\phi_{ex}$, the approximate solution $\\phi_h$ and the error function $|\\phi_{ex} - \\phi_h|$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from sympy import lambdify\n", + "\n", + "from utils import plot\n", + "\n", + "phi_ex_fun = lambdify(domain.coordinates, phi_ex)\n", + "error = lambda x, y: abs(phi_ex_fun(x, y) - phi_h_FemField(x, y))\n", + "\n", + "plot(gridsize_x = 100, \n", + " gridsize_y = 100, \n", + " title = r'Approximation of Solution $\\phi$', \n", + " funs = [phi_ex_fun, phi_h_FemField, error], \n", + " titles = [r'$\\phi_{ex}(x,y)$', r'$\\phi_h(x,y)$', r'$|(\\phi_{ex}-\\phi_h)(x,y)|$'],\n", + " surface_plot = True\n", + ")" + ] + } + ], + "metadata": {}, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/exercises/time_harmonic_maxwell_2d_I.ipynb b/exercises/time_harmonic_maxwell_2d_I.ipynb new file mode 100644 index 0000000..98be3e1 --- /dev/null +++ b/exercises/time_harmonic_maxwell_2d_I.ipynb @@ -0,0 +1,291 @@ +{ + "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 time-harmonic Maxwell equation on the Unit Square\n", + "## ... using Psydac's bilinear form interface\n", + "\n", + "In this exercise we write a solver for the 2D time-harmonic Maxwell equation on the unit square using Psydac's bilinear form interface.\n", + "\n", + "\\begin{align*}\n", + " \\boldsymbol{curl}curl\\boldsymbol{E} - \\omega^2\\boldsymbol{E} &= \\boldsymbol{F} \\quad \\text{ in }\\Omega=]0,1[^2, \\\\\n", + " \\boldsymbol{n}\\times\\boldsymbol{E} &= 0 \\quad \\text{ on }\\partial\\Omega,\n", + "\\end{align*}\n", + "for the specific choice\n", + "\\begin{align*}\n", + " \\omega &= 1.5 \\\\\n", + " \\alpha &= -\\omega^2 \\\\\n", + " \\boldsymbol{F}(x, y) &= \\left(\\begin{matrix}\n", + " (\\alpha+\\pi^2)\\sin(\\pi y) - \\pi^2\\cos(\\pi x)\\sin(\\pi y) \\\\\n", + " (\\alpha+\\pi^2)\\sin(\\pi x)\\cos(\\pi y)\n", + " \\end{matrix}\\right)\n", + "\\end{align*}\n", + "\n", + "## The Variational Formulation\n", + "\n", + "The corresponding variational formulation reads\n", + "\n", + "\\begin{align*}\n", + " \\text{Find }\\boldsymbol{E}\\in V\\coloneqq H_0(curl;\\Omega)\\text{ such that } \\\\\n", + " a(\\boldsymbol{E}, \\boldsymbol{v}) = l(\\boldsymbol{v})\\quad\\forall\\ \\boldsymbol{v}\\in V\n", + "\\end{align*}\n", + "\n", + "where \n", + "- $a(\\boldsymbol{E},\\boldsymbol{v}) \\coloneqq \\int_{\\Omega} (\\nabla\\times\\boldsymbol{E})(\\nabla\\times\\boldsymbol{v}) + \\alpha \\boldsymbol{E}\\cdot\\boldsymbol{v} ~ d\\Omega$ ,\n", + "- $l(\\boldsymbol{v}) := \\int_{\\Omega} \\boldsymbol{F}\\cdot\\boldsymbol{v} ~ d\\Omega$." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Formal Model" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from sympde.calculus import dot, curl\n", + "from sympde.topology import VectorFunctionSpace, elements_of, Square\n", + "from sympde.expr.expr import LinearForm, BilinearForm, integral\n", + "\n", + "domain = Square('S', bounds1=(0, 1), bounds2=(0, 1))\n", + "\n", + "V = VectorFunctionSpace('V', domain, kind='hcurl')\n", + "\n", + "E, v = elements_of(V, names='E, v')\n", + "\n", + "# bilinear form\n", + "omega = 1.5\n", + "alpha = -omega**2\n", + "expr1 = dot(E, v)\n", + "a = BilinearForm((E, v), integral(domain, expr1))\n", + "\n", + "# linear form\n", + "from sympy import pi, sin, cos, Tuple, Matrix\n", + "x,y = domain.coordinates\n", + "F = Tuple( (alpha + pi**2) * sin(pi*y) - pi**2 * sin(pi*y) * cos(pi*x),\n", + " (alpha + pi**2) * sin(pi*x) * cos(pi*y) )\n", + "\n", + "expr2 = dot(F, v)\n", + "l = LinearForm(v, integral(domain, expr2))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Discretization\n", + "\n", + "We will use Psydac to discretize the problem." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from psydac.api.discretization import discretize\n", + "from psydac.api.settings import PSYDAC_BACKEND_GPYCCEL\n", + "\n", + "backend = PSYDAC_BACKEND_GPYCCEL" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ncells = [16, 16] # Bspline cells\n", + "degree = [3, 3] # Bspline degree" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "domain_h = discretize(domain, ncells=ncells, periodic=[False, False])\n", + "V_h = discretize(V, domain_h, degree=degree)\n", + "\n", + "a_h = discretize(a, domain_h, (V_h, V_h), backend=backend)\n", + "l_h = discretize(l, domain_h, V_h, backend=backend)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Boundary Conditions\n", + "\n", + "We choose to apply a projection method. For that matter, we construct the projection matrix $\\mathbb{P}_0$ and its counterpart $\\mathbb{P}_{\\Gamma} = \\mathbb{I} - \\mathbb{P}_0$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from utils import HcurlBoundaryProjector2D\n", + "from psydac.linalg.basic import IdentityOperator\n", + "\n", + "P_Hcurl = HcurlBoundaryProjector2D(V, V_h.coeff_space)\n", + "\n", + "I1 = IdentityOperator(V_h.coeff_space)\n", + "P_Gamma = I1 - P_Hcurl" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "A = a_h.assemble()\n", + "f = l_h.assemble()\n", + "\n", + "A_bc = P_Hcurl @ A @ P_Hcurl + P_Gamma\n", + "f_bc = P_Hcurl @ f" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Solving the PDE" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import time\n", + "\n", + "from psydac.linalg.solvers import inverse\n", + "\n", + "tol = 1e-9\n", + "maxiter = 1000\n", + "\n", + "A_bc_inv = inverse(A_bc, 'cg', tol=tol, maxiter=maxiter)\n", + "\n", + "t0 = time.time()\n", + "E_h = A_bc_inv @ f_bc\n", + "t1 = time.time()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Computing the error norm\n", + "\n", + "As the analytical solution is available, we want to compute the $L^2$ norm of the error.\n", + "In this example, the analytical solution is given by\n", + "\n", + "$$\n", + "\\boldsymbol{E}_{ex}(x, y) = \\left(\\begin{matrix} \n", + " \\sin(\\pi y) \\\\\n", + " \\sin(\\pi x)\\cos(\\pi y)\n", + " \\end{matrix}\\right)\n", + "$$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from psydac.fem.basic import FemField\n", + "\n", + "from sympde.expr import Norm\n", + "\n", + "E_ex = Tuple(sin(pi*y), sin(pi*x)*cos(pi*y))\n", + "E_h_FemField = FemField(V_h, E_h)\n", + "\n", + "error = Matrix([E[0] - E_ex[0], E[1] - E_ex[1]])\n", + "\n", + "# create the formal Norm object\n", + "l2norm = Norm(error, domain, kind='l2')\n", + "\n", + "# discretize the norm\n", + "l2norm_h = discretize(l2norm, domain_h, V_h, backend=backend)\n", + "\n", + "# assemble the norm\n", + "l2_error = l2norm_h.assemble(E=E_h_FemField)\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( '> CG info :: ',A_bc_inv.get_info() )\n", + "print( '> L2 error :: {:.2e}'.format( l2_error ) )\n", + "print( '' )\n", + "print( '> Solution time :: {:.3g}'.format( t1-t0 ) )" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Visualization\n", + "\n", + "We plot the true solution $\\boldsymbol{E}_{ex}$, the approximate solution $\\boldsymbol{E}_h$ and the error function $|\\boldsymbol{E}_{ex} - \\boldsymbol{E}_h|$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from sympy import lambdify\n", + "\n", + "from utils import plot\n", + "\n", + "E_ex_x = lambdify((x, y), E_ex[0])\n", + "E_ex_y = lambdify((x, y), E_ex[1])\n", + "E_h_x = E_h_FemField[0]\n", + "E_h_y = E_h_FemField[1]\n", + "error_x = lambda x, y: abs(E_ex_x(x, y) - E_h_x(x, y))\n", + "error_y = lambda x, y: abs(E_ex_y(x, y) - E_h_y(x, y))\n", + "\n", + "plot(gridsize_x = 100, \n", + " gridsize_y = 100, \n", + " title = r'Approximation of Solution $\\boldsymbol{E}$, first component', \n", + " funs = [E_ex_x, E_h_x, error_x], \n", + " titles = [r'$(\\boldsymbol{E}_{ex})_1(x,y)$', r'$(\\boldsymbol{E}_h)_1(x,y)$', r'$|(\\boldsymbol{E}_{ex}-\\boldsymbol{E}_h)_1(x,y)|$'],\n", + " surface_plot = True\n", + ")\n", + "\n", + "plot(gridsize_x = 100,\n", + " gridsize_y = 100,\n", + " title = r'Approximation of Solution $\\boldsymbol{E}$, second component',\n", + " funs = [E_ex_y, E_h_y, error_y],\n", + " titles = [r'$(\\boldsymbol{E}_{ex})_2(x,y)$', r'$(\\boldsymbol{E}_h)_2(x,y)$', r'$|(\\boldsymbol{E}_{ex}-\\boldsymbol{E}_h)_2(x,y)|$'],\n", + " surface_plot = True\n", + ")" + ] + } + ], + "metadata": {}, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/exercises/time_harmonic_maxwell_2d_II.ipynb b/exercises/time_harmonic_maxwell_2d_II.ipynb new file mode 100644 index 0000000..abf33ef --- /dev/null +++ b/exercises/time_harmonic_maxwell_2d_II.ipynb @@ -0,0 +1,299 @@ +{ + "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 time-harmonic Maxwell equation on the Unit Square\n", + "## ... using Psydac's de Rham interface\n", + "\n", + "In this exercise we write a solver for the 2D time-harmonic Maxwell equation on the unit square using Psydac's de Rham interface.\n", + "\n", + "\\begin{align*}\n", + " \\boldsymbol{curl}curl\\boldsymbol{E} - \\omega^2\\boldsymbol{E} &= \\boldsymbol{F} \\quad \\text{ in }\\Omega=]0,1[^2, \\\\\n", + " \\boldsymbol{n}\\times\\boldsymbol{E} &= 0 \\quad \\text{ on }\\partial\\Omega,\n", + "\\end{align*}\n", + "for the specific choice\n", + "\\begin{align*}\n", + " \\omega &= 1.5 \\\\\n", + " \\alpha &= -\\omega^2 \\\\\n", + " \\boldsymbol{F}(x, y) &= \\left(\\begin{matrix}\n", + " (\\alpha+\\pi^2)\\sin(\\pi y) - \\pi^2\\cos(\\pi x)\\sin(\\pi y) \\\\\n", + " (\\alpha+\\pi^2)\\sin(\\pi x)\\cos(\\pi y)\n", + " \\end{matrix}\\right)\n", + "\\end{align*}\n", + "\n", + "## The Variational Formulation\n", + "\n", + "The corresponding variational formulation reads\n", + "\n", + "\\begin{align*}\n", + " \\text{Find }\\boldsymbol{E}\\in V\\coloneqq H_0(curl;\\Omega)\\text{ such that } \\\\\n", + " a(\\boldsymbol{E}, \\boldsymbol{v}) = l(\\boldsymbol{v})\\quad\\forall\\ \\boldsymbol{v}\\in V\n", + "\\end{align*}\n", + "\n", + "where \n", + "- $a(\\boldsymbol{E},\\boldsymbol{v}) \\coloneqq \\int_{\\Omega} (\\nabla\\times\\boldsymbol{E})(\\nabla\\times\\boldsymbol{v}) + \\alpha \\boldsymbol{E}\\cdot\\boldsymbol{v} ~ d\\Omega$ ,\n", + "- $l(\\boldsymbol{v}) := \\int_{\\Omega} \\boldsymbol{F}\\cdot\\boldsymbol{v} ~ d\\Omega$." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Discrete Model using de Rham objects" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from sympde.calculus import dot\n", + "from sympde.topology import elements_of, Square, Derham\n", + "from sympde.expr.expr import LinearForm, BilinearForm, integral\n", + "\n", + "domain = Square('S', bounds1=(0, 1), bounds2=(0, 1))\n", + "derham = Derham(domain, sequence=['h1', 'hcurl', 'l2'])\n", + "\n", + "V1 = derham.V1\n", + "V2 = derham.V2\n", + "\n", + "u1, v1 = elements_of(V1, names='u1, v1')\n", + "u2, v2 = elements_of(V2, names='u2, v2')\n", + "\n", + "# bilinear forms corresponding to V1 and V2 mass matrices\n", + "m1 = BilinearForm((u1, v1), integral(domain, dot(u1, v1)))\n", + "m2 = BilinearForm((u2, v2), integral(domain, u2 * v2))\n", + "\n", + "# linear form corresponding to the rhs\n", + "from sympy import pi, sin, cos, Tuple, Matrix\n", + "omega = 1.5\n", + "alpha = -omega**2\n", + "x,y = domain.coordinates\n", + "F = Tuple( (alpha + pi**2) * sin(pi*y) - pi**2 * sin(pi*y) * cos(pi*x),\n", + " (alpha + pi**2) * sin(pi*x) * cos(pi*y) )\n", + "\n", + "expr = dot(F, v1)\n", + "l = LinearForm(v1, integral(domain, expr))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from psydac.api.discretization import discretize\n", + "from psydac.api.settings import PSYDAC_BACKEND_GPYCCEL\n", + "\n", + "backend = PSYDAC_BACKEND_GPYCCEL" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ncells = [16, 16] # Bspline cells\n", + "degree = [3, 3] # Bspline degree" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "domain_h = discretize(domain, ncells=ncells, periodic=[False, False])\n", + "derham_h = discretize(derham, domain_h, degree=degree)\n", + "\n", + "V1_h = derham_h.V1\n", + "V2_h = derham_h.V2\n", + "\n", + "# Exterior Derivative operators (grad and curl)\n", + "G, C = derham_h.derivatives_as_matrices\n", + "\n", + "# Mass matrices\n", + "m1_h = discretize(m1, domain_h, (V1_h, V1_h), backend=backend)\n", + "m2_h = discretize(m2, domain_h, (V2_h, V2_h), backend=backend)\n", + "\n", + "M1 = m1_h.assemble()\n", + "M2 = m2_h.assemble()\n", + "\n", + "# System Matrix A\n", + "A = M1\n", + "\n", + "# rhs vector f\n", + "l_h = discretize(l, domain_h, V1_h, backend=backend)\n", + "f = l_h.assemble()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Boundary Conditions\n", + "\n", + "We choose to apply a projection method. For that matter, we construct the projection matrix $\\mathbb{P}_0$ and its counterpart $\\mathbb{P}_{\\Gamma} = \\mathbb{I} - \\mathbb{P}_0$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from utils import HcurlBoundaryProjector2D\n", + "from psydac.linalg.basic import IdentityOperator\n", + "\n", + "P_Hcurl = HcurlBoundaryProjector2D(V1, V1_h.coeff_space)\n", + "\n", + "I1 = IdentityOperator(V1_h.coeff_space)\n", + "P_Gamma = I1 - P_Hcurl" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "A_bc = P_Hcurl @ A @ P_Hcurl + P_Gamma\n", + "f_bc = P_Hcurl @ f" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Solving the PDE" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import time\n", + "\n", + "from psydac.linalg.solvers import inverse\n", + "\n", + "tol = 1e-9\n", + "maxiter = 1000\n", + "\n", + "A_bc_inv = inverse(A_bc, 'cg', tol=tol, maxiter=maxiter)\n", + "\n", + "t0 = time.time()\n", + "E_h = A_bc_inv @ f_bc\n", + "t1 = time.time()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Computing the error norm\n", + "\n", + "As the analytical solution is available, we want to compute the $L^2$ norm of the error.\n", + "In this example, the analytical solution is given by\n", + "\n", + "$$\n", + "\\boldsymbol{E}_{ex}(x, y) = \\left(\\begin{matrix} \n", + " \\sin(\\pi y) \\\\\n", + " \\sin(\\pi x)\\cos(\\pi y)\n", + " \\end{matrix}\\right)\n", + "$$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from psydac.fem.basic import FemField\n", + "\n", + "from sympde.expr import Norm\n", + "\n", + "E_ex = Tuple(sin(pi*y), sin(pi*x)*cos(pi*y))\n", + "E_h_FemField = FemField(V1_h, E_h)\n", + "\n", + "error = Matrix([u1[0] - E_ex[0], u1[1] - E_ex[1]])\n", + "\n", + "# create the formal Norm object\n", + "l2norm = Norm(error, domain, kind='l2')\n", + "\n", + "# discretize the norm\n", + "l2norm_h = discretize(l2norm, domain_h, V1_h, backend=backend)\n", + "\n", + "# assemble the norm\n", + "l2_error = l2norm_h.assemble(u1=E_h_FemField)\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( '> CG info :: ',A_bc_inv.get_info() )\n", + "print( '> L2 error :: {:.2e}'.format( l2_error ) )\n", + "print( '' )\n", + "print( '> Solution time :: {:.3g}'.format( t1-t0 ) )" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Visualization\n", + "\n", + "We plot the true solution $\\boldsymbol{E}_{ex}$, the approximate solution $\\boldsymbol{E}_h$ and the error function $|\\boldsymbol{E}_{ex} - \\boldsymbol{E}_h|$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from sympy import lambdify\n", + "\n", + "from utils import plot\n", + "\n", + "E_ex_x = lambdify((x, y), E_ex[0])\n", + "E_ex_y = lambdify((x, y), E_ex[1])\n", + "E_h_x = E_h_FemField[0]\n", + "E_h_y = E_h_FemField[1]\n", + "error_x = lambda x, y: abs(E_ex_x(x, y) - E_h_x(x, y))\n", + "error_y = lambda x, y: abs(E_ex_y(x, y) - E_h_y(x, y))\n", + "\n", + "plot(gridsize_x = 100, \n", + " gridsize_y = 100, \n", + " title = r'Approximation of Solution $\\boldsymbol{E}$, first component', \n", + " funs = [E_ex_x, E_h_x, error_x], \n", + " titles = [r'$(\\boldsymbol{E}_{ex})_1(x,y)$', r'$(\\boldsymbol{E}_h)_1(x,y)$', r'$|(\\boldsymbol{E}_{ex}-\\boldsymbol{E}_h)_1(x,y)|$'],\n", + " surface_plot = True\n", + ")\n", + "\n", + "plot(gridsize_x = 100,\n", + " gridsize_y = 100,\n", + " title = r'Approximation of Solution $\\boldsymbol{E}$, second component',\n", + " funs = [E_ex_y, E_h_y, error_y],\n", + " titles = [r'$(\\boldsymbol{E}_{ex})_2(x,y)$', r'$(\\boldsymbol{E}_h)_2(x,y)$', r'$|(\\boldsymbol{E}_{ex}-\\boldsymbol{E}_h)_2(x,y)|$'],\n", + " surface_plot = True\n", + ")" + ] + } + ], + "metadata": {}, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/exercises/utils.py b/exercises/utils.py new file mode 100644 index 0000000..27abb27 --- /dev/null +++ b/exercises/utils.py @@ -0,0 +1,382 @@ +import numpy as np +import matplotlib.pyplot as plt +from matplotlib import cm, colors + +from scipy.sparse import dia_matrix +from scipy.sparse.linalg import eigsh + +from sympde.expr import EssentialBC, BilinearForm, integral +from sympde.topology import element_of, elements_of, Line, Derham +from sympde.topology.datatype import H1Space, HcurlSpace, L2Space + +from psydac.api.discretization import discretize +from psydac.api.essential_bc import apply_essential_bc +from psydac.linalg.basic import LinearOperator, Vector +from psydac.linalg.block import BlockVectorSpace, BlockLinearOperator +from psydac.linalg.direct_solvers import BandedSolver +from psydac.linalg.kron import KroneckerLinearSolver +from psydac.linalg.stencil import StencilVectorSpace + + +def get_unit_vector(v, n1, n2, n3, pads1, pads2, pads3): + + v *= 0.0 + if n3 is None: + assert pads3 is None + if n2 is None: + raise NotImplementedError('This get_unit_vector method is only implemented in 2D.') + else: + v._data[pads1+n1, pads2+n2] = 1. + else: + raise NotImplementedError('This get_unit_vector method is only implemented in 2D.') + + return v + +def toarray(A): + """Obtain a numpy array representation of a LinearOperator (which has not implemented toarray()).""" + assert isinstance(A, LinearOperator) + + W = A.codomain + + W_is_block = True if isinstance(W, BlockVectorSpace) else False + if not W_is_block: + assert isinstance(W, StencilVectorSpace) + + A_arr = np.zeros(A.shape, dtype="float64") + w = W.zeros() + At = A.T + + if W_is_block: + codomain_blocks = [W[i] for i in range(W.n_blocks)] + else: + codomain_blocks = (W, ) + + start_index = 0 + for k, Wk in enumerate(codomain_blocks): + w *= 0. + v = Wk.zeros() + if len(Wk.npts) == 2: + npts1, npts2 = Wk.npts + pads1, pads2 = Wk.pads + for n1 in range(npts1): + for n2 in range(npts2): + e_n1_n2 = get_unit_vector(v, n1, n2, None, pads1, pads2, None) + if W_is_block: + w[k] = e_n1_n2 + e_n1_n2 = w + A_n1_n2 = At @ e_n1_n2 + A_arr[start_index + n1*npts2+n2, :] = A_n1_n2.toarray() + else: + raise NotImplementedError('This toarray method is currently only implemented in 2D.') + start_index += Wk.dimension + + return A_arr + +class H1BoundaryProjector2D(LinearOperator): + def __init__(self, V0, V0_vs, periodic=[False, False]): + + assert all([isinstance(P, bool) for P in periodic]) + + self._domain = V0_vs + self._codomain = V0_vs + self._space = V0 + self._periodic = periodic + + self._BC = self._get_BC() + + def copy(self): + return H1BoundaryProjector2D(self._space, self._domain, self._periodic) + + def _get_BC(self): + + periodic = self._periodic + if all([P == True for P in periodic]): + return None + + space = self._space + u = element_of(space, name='u') + bcs = [EssentialBC(u, 0, side, position=0) for side in space.domain.boundary] + + bcs_x = [bcs[0], bcs[1]] if periodic[0] == False else [] + bcs_y = [bcs[2], bcs[3]] if periodic[1] == False else [] + + BC = bcs_x + bcs_y + + return BC + + @property + def domain(self): + return self._domain + + @property + def codomain(self): + return self._codomain + + @property + def dtype(self): + return None + + def tosparse(self): + raise NotImplementedError + + def toarray(self): + return toarray(self) + + def transpose(self, conjugate=False): + return self + + def dot(self, v, out=None): + BC = self._BC + if out is not None: + assert isinstance(out, Vector) + assert out.space is self.codomain + else: + out = self.codomain.zeros() + + v.copy(out=out) + apply_essential_bc(out, *BC) + return out + +class HcurlBoundaryProjector2D(LinearOperator): + def __init__(self, V1, V1_vs, periodic=[False, False]): + + assert all([isinstance(P, bool) for P in periodic]) + + self._domain = V1_vs + self._codomain = V1_vs + self._space = V1 + self._periodic = periodic + + self._BC = self._get_BC() + + def copy(self): + return HcurlBoundaryProjector2D(self._space, self._domain, self._periodic) + + def _get_BC(self): + + periodic = self._periodic + if all([P == True for P in periodic]): + return None + + space = self._space + u = element_of(space, name='u') + bcs = [EssentialBC(u, 0, side, position=0) for side in space.domain.boundary] + + bcs_x = [bcs[0], bcs[1]] if periodic[0] == False else [] + bcs_y = [bcs[2], bcs[3]] if periodic[1] == False else [] + + BC = bcs_x + bcs_y + + bcs_x = [bcs[2], bcs[3]] if periodic[1] == False else [] + bcs_y = [bcs[0], bcs[1]] if periodic[0] == False else [] + + BC = [bcs_x, bcs_y] + + return BC + + @property + def domain(self): + return self._domain + + @property + def codomain(self): + return self._codomain + + @property + def dtype(self): + return None + + def tosparse(self): + raise NotImplementedError + + def toarray(self): + return toarray(self) + + def transpose(self, conjugate=False): + return self + + def dot(self, v, out=None): + BC = self._BC + if out is not None: + assert isinstance(out, Vector) + assert out.space is self.codomain + else: + out = self.codomain.zeros() + + v.copy(out=out) + for outi, BCi in zip(out, BC): + apply_essential_bc(outi, *BCi) + return out + +def plot(gridsize_x, gridsize_y, title, funs, titles, surface_plot=False): + + x = np.linspace(0, 1, gridsize_x+1) + y = np.linspace(0, 1, gridsize_y+1) + xx, yy = np.meshgrid(x, y) + vals = [[] for _ in funs] + for i, fun in enumerate(funs): + for xi, yi in zip(xx, yy): + vals[i].append([fun(xii, yii) for xii, yii in zip(xi, yi)]) + vals[i] = np.array(vals[i]) + + n_plots = len(funs) + if n_plots > 1: + assert n_plots == len(titles) + else: + if titles is not None: + print('Warning [plot]: will discard argument titles for a single plot') + + fig = plt.figure(figsize=(2.6+4.8*n_plots, 4.8)) + fig.suptitle(title, fontsize=14) + + for i in range(n_plots): + vmin = np.min(vals[i]) + vmax = np.max(vals[i]) + cnorm = colors.Normalize(vmin=vmin, vmax=vmax) + ax = fig.add_subplot(1, n_plots, i+1) + ax.contourf(xx, yy, vals[i], 50, norm=cnorm, cmap='viridis') + ax.axis('equal') + fig.colorbar(cm.ScalarMappable(norm=cnorm, cmap='viridis'), ax=ax, pad=0.05) + ax.set_xlabel( r'$x$', rotation='horizontal' ) + ax.set_ylabel( r'$y$', rotation='horizontal' ) + if n_plots > 1: + ax.set_title ( titles[i] ) + plt.show() + + if surface_plot: + fig = plt.figure(figsize=(2.6+4.8*n_plots, 4.8)) + fig.suptitle(title+' -- surface', fontsize=14) + + for i in range(n_plots): + vmin = np.min(vals[i]) + vmax = np.max(vals[i]) + cnorm = colors.Normalize(vmin=vmin, vmax=vmax) + ax = fig.add_subplot(1, n_plots, i+1, projection='3d') + ax.plot_surface(xx, yy, vals[i], norm=cnorm, cmap='viridis', + linewidth=0, antialiased=False) + fig.colorbar(cm.ScalarMappable(norm=cnorm, cmap='viridis'), ax=ax, pad=0.05) + ax.set_xlabel( r'$x$', rotation='horizontal' ) + ax.set_ylabel( r'$y$', rotation='horizontal' ) + if n_plots > 1: + ax.set_title ( titles[i] ) + plt.show() + +def to_bnd(A): + + dmat = dia_matrix(A.toarray(), dtype=A.dtype) + la = abs(dmat.offsets.min()) + ua = dmat.offsets.max() + cmat = dmat.tocsr() + + A_bnd = np.zeros((1+ua+2*la, cmat.shape[1]), A.dtype) + + for i,j in zip(*cmat.nonzero()): + A_bnd[la+ua+i-j, j] = cmat[i,j] + + return A_bnd, la, ua + +def matrix_to_bandsolver(A): + A.remove_spurious_entries() + A_bnd, la, ua = to_bnd(A) + return BandedSolver(ua, la, A_bnd) + +def get_M1_block_kron_solver_2D(V1, ncells, degree, periodic): + """ + Given a 2D DeRham sequenece (V0 = H(grad) --grad--> V1 = H(curl) --curl--> V2 = L2) + discreticed using ncells, degree and periodic, + + domain = Square('C', bounds1=(0, 1), bounds2=(0, 1)) + derham = Derham(domain) + domain_h = discretize(domain, ncells=ncells, periodic=periodic, comm=comm) + derham_h = discretize(derham, domain_h, degree=degree), + + returns the inverse of the mass matrix M1 as a BlockLinearOperator consisting of two KroneckerLinearSolvers on the diagonal. + """ + # assert 3D + assert len(ncells) == 2 + assert len(degree) == 2 + assert len(periodic) == 2 + + # 1D domain to be discreticed using the respective values of ncells, degree, periodic + domain_1d = Line('L', bounds=(0,1)) + derham_1d = Derham(domain_1d) + + # storage for the 1D mass matrices + M0_matrices = [] + M1_matrices = [] + + # assembly of the 1D mass matrices + for (n, p, P) in zip(ncells, degree, periodic): + + domain_1d_h = discretize(domain_1d, ncells=[n], periodic=[P]) + derham_1d_h = discretize(derham_1d, domain_1d_h, degree=[p]) + + u_1d_0, v_1d_0 = elements_of(derham_1d.V0, names='u_1d_0, v_1d_0') + u_1d_1, v_1d_1 = elements_of(derham_1d.V1, names='u_1d_1, v_1d_1') + + a_1d_0 = BilinearForm((u_1d_0, v_1d_0), integral(domain_1d, u_1d_0 * v_1d_0)) + a_1d_1 = BilinearForm((u_1d_1, v_1d_1), integral(domain_1d, u_1d_1 * v_1d_1)) + + a_1d_0_h = discretize(a_1d_0, domain_1d_h, (derham_1d_h.V0, derham_1d_h.V0)) + a_1d_1_h = discretize(a_1d_1, domain_1d_h, (derham_1d_h.V1, derham_1d_h.V1)) + + M_1d_0 = a_1d_0_h.assemble() + M_1d_1 = a_1d_1_h.assemble() + + M0_matrices.append(M_1d_0) + M1_matrices.append(M_1d_1) + + V1_1 = V1[0] + V1_2 = V1[1] + + B1_mat = [M1_matrices[0], M0_matrices[1]] + B2_mat = [M0_matrices[0], M1_matrices[1]] + + B1_solvers = [matrix_to_bandsolver(Ai) for Ai in B1_mat] + B2_solvers = [matrix_to_bandsolver(Ai) for Ai in B2_mat] + + B1_kron_inv = KroneckerLinearSolver(V1_1, V1_1, B1_solvers) + B2_kron_inv = KroneckerLinearSolver(V1_2, V1_2, B2_solvers) + + M1_block_kron_solver = BlockLinearOperator(V1, V1, ((B1_kron_inv, None), + (None, B2_kron_inv))) + + return M1_block_kron_solver + +def get_eigenvalues(nb_eigs, sigma, A_m, M_m): + """ + Compute the eigenvalues of the matrix A close to sigma and right-hand-side M + Function seen and adapted from >>> psydac_dev/psydac/feec/multipatch/examples/hcurl_eigen_pbms_conga_2d.py <<< + (Commit a748a4d8c1569a8765f6688d228f65ea6073c252) + + Parameters + ---------- + nb_eigs : int + Number of eigenvalues to compute + sigma : float + Value close to which the eigenvalues are computed + A_m : sparse matrix + Matrix A + M_m : sparse matrix + Matrix M + """ + + print() + print('----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ') + print( + 'computing {0} eigenvalues (and eigenvectors) close to sigma={1} with scipy.sparse.eigsh...'.format(nb_eigs, sigma)) + mode = 'normal' + which = 'LM' + ncv = 4 * nb_eigs + max_shape_splu = 24000 + if A_m.shape[0] >= max_shape_splu: + raise ValueError(f'Matrix too large.') + + eigenvalues, eigenvectors = eigsh( + A_m, k=nb_eigs, M=M_m, sigma=sigma, mode=mode, which=which, ncv=ncv) + + print("done: eigenvalues found: " + repr(eigenvalues)) + print('----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- -----') + print() + + return eigenvalues, eigenvectors