From 29cc8e153f6149501dfa913c39ac64855b4979dd Mon Sep 17 00:00:00 2001 From: Martin Campos Pinto Date: Tue, 11 Mar 2025 15:34:17 +0100 Subject: [PATCH 01/16] new exercise files --- _toc.yml | 4 + exercises/harmonic-feec.ipynb | 274 ++++++++++++++++++++++++++++++++++ exercises/harmonic-fem.ipynb | 274 ++++++++++++++++++++++++++++++++++ 3 files changed, 552 insertions(+) create mode 100644 exercises/harmonic-feec.ipynb create mode 100644 exercises/harmonic-fem.ipynb diff --git a/_toc.yml b/_toc.yml index 1db5e02..cd82ca3 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/harmonic-feec.ipynb b/exercises/harmonic-feec.ipynb new file mode 100644 index 0000000..65af55d --- /dev/null +++ b/exercises/harmonic-feec.ipynb @@ -0,0 +1,274 @@ +{ + "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": { + "language_info": { + "name": "python" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/exercises/harmonic-fem.ipynb b/exercises/harmonic-fem.ipynb new file mode 100644 index 0000000..7d3f6d9 --- /dev/null +++ b/exercises/harmonic-fem.ipynb @@ -0,0 +1,274 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "9f28f9af", + "metadata": {}, + "source": [ + "# Harmonic fields with H1 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": { + "language_info": { + "name": "python" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} From 186494b794237a239d5f1405e1d3bbaaf5fe0e53 Mon Sep 17 00:00:00 2001 From: Martin Campos Pinto Date: Tue, 11 Mar 2025 16:25:09 +0100 Subject: [PATCH 02/16] pbm formulation --- .gitignore | 2 ++ exercises/harmonic-feec.ipynb | 6 +----- exercises/harmonic-fem.ipynb | 30 ++++++++++++++++++------------ 3 files changed, 21 insertions(+), 17 deletions(-) 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/exercises/harmonic-feec.ipynb b/exercises/harmonic-feec.ipynb index 65af55d..d7da6a7 100644 --- a/exercises/harmonic-feec.ipynb +++ b/exercises/harmonic-feec.ipynb @@ -264,11 +264,7 @@ ] } ], - "metadata": { - "language_info": { - "name": "python" - } - }, + "metadata": {}, "nbformat": 4, "nbformat_minor": 5 } diff --git a/exercises/harmonic-fem.ipynb b/exercises/harmonic-fem.ipynb index 7d3f6d9..02a87d9 100644 --- a/exercises/harmonic-fem.ipynb +++ b/exercises/harmonic-fem.ipynb @@ -7,12 +7,13 @@ "source": [ "# Harmonic fields with H1 fem spaces\n", "\n", - "In this example we .... consider the vector Poisson equation with homogeneous Dirichlet boundary conditions:\n", + "In this exercise we consider the following vector problem with homogeneous boundary conditions in H(curl):\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", + " \\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", @@ -23,15 +24,24 @@ "$$\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", + " 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}_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$." + "- $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" ] }, { @@ -264,11 +274,7 @@ ] } ], - "metadata": { - "language_info": { - "name": "python" - } - }, + "metadata": {}, "nbformat": 4, "nbformat_minor": 5 } From 7fcf870fc49e545d309b72e0c7be4db80d976a1b Mon Sep 17 00:00:00 2001 From: jowezarek Date: Mon, 24 Mar 2025 15:35:05 +0100 Subject: [PATCH 03/16] add sheet1 exercise 2.2 --- exercises/poisson_2d_I.ipynb | 311 +++++++++++++++++++++++++++++++++++ exercises/utils.py | 197 ++++++++++++++++++++++ 2 files changed, 508 insertions(+) create mode 100644 exercises/poisson_2d_I.ipynb create mode 100644 exercises/utils.py diff --git a/exercises/poisson_2d_I.ipynb b/exercises/poisson_2d_I.ipynb new file mode 100644 index 0000000..0c2ec2d --- /dev/null +++ b/exercises/poisson_2d_I.ipynb @@ -0,0 +1,311 @@ +{ + "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", + "\n", + "In the 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", + "\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\n", + "from sympde.topology import ScalarFunctionSpace\n", + "\n", + "from sympy import pi, sin\n", + "\n", + "domain = Square('S', bounds1=(0, 1), bounds2=(0, 1))\n", + "\n", + "V = ScalarFunctionSpace('V', domain)\n", + "x,y = domain.coordinates\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", + "f = 8 * (pi**2) * sin(2 * pi * x) * sin(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" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ncells = [8, 8] # Bspline cells\n", + "degree = [2, 2] # Bspline degree" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from psydac.api.settings import PSYDAC_BACKEND_GPYCCEL\n", + "\n", + "backend = PSYDAC_BACKEND_GPYCCEL\n", + "\n", + "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_2\n", + "from psydac.linalg.basic import IdentityOperator\n", + "\n", + "P_0 = H1BoundaryProjector2D_2(V, V_h.vector_space)\n", + "\n", + "I_0 = IdentityOperator(V_h.vector_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": [ + "from psydac.linalg.solvers import inverse\n", + "\n", + "tol = 1e-12\n", + "maxiter = 1000\n", + "\n", + "phi_h = inverse(A_bc, 'cg', tol=tol, maxiter=maxiter) @ f_bc" + ] + }, + { + "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) = \\sin(2\\pi x)\\sin(2\\pi y)\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Computing the $L^2$ norm" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from sympde.expr import find, EssentialBC, Norm, SemiNorm\n", + "from psydac.fem.basic import FemField\n", + "\n", + "phi_ex = sin(2*pi*x) * sin(2*pi*y)\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)\n", + "\n", + "# assemble the norm\n", + "l2_error = l2norm_h.assemble(phi=FemField(V_h, phi_h))\n", + "\n", + "# print the result\n", + "print(l2_error)" + ] + }, + { + "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", + "h1_error = h1norm_h.assemble(phi=FemField(V_h, phi_h))\n", + "\n", + "# print the result\n", + "print(h1_error)" + ] + }, + { + "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=FemField(V_h, phi_h))\n", + "\n", + "# print the result\n", + "print(h1_error)" + ] + }, + { + "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 utils import plot\n", + "from sympy import lambdify\n", + "\n", + "phi_ex_fun = lambdify(domain.coordinates, phi_ex)\n", + "phi_h_fun = FemField(V_h, phi_h)\n", + "error = lambda x, y: abs(phi_ex_fun(x, y) - phi_h_fun(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_fun, 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": 2 +} diff --git a/exercises/utils.py b/exercises/utils.py new file mode 100644 index 0000000..c777105 --- /dev/null +++ b/exercises/utils.py @@ -0,0 +1,197 @@ +import h5py +import numpy as np +import matplotlib.pyplot as plt +from matplotlib import cm, colors + +from sympde.expr import EssentialBC +from sympde.topology import element_of + +from psydac.api.essential_bc import apply_essential_bc +from psydac.linalg.basic import LinearOperator, Vector +from psydac.linalg.block import BlockVectorSpace +from psydac.linalg.stencil import StencilVectorSpace + + +def get_unit_vector(v, n1, n2, pads1, pads2): + + v *= 0.0 + v._data[pads1+n1, pads2+n2] = 1. + + 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() + 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, pads1, pads2) + 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() + start_index += Wk.dimension + + return A_arr + +def plot(gridsize_x, title, funs, titles, gridsize_y=None, surface_plot=False): + + if gridsize_y is None: + x = np.linspace(0, 1, gridsize_x+1) + vals = [[] for fun in funs] + for i, fun in enumerate(funs): + for xi in x: + vals[i].append(fun(xi)) + + n_plots = len(funs) + if n_plots > 1: + assert n_plots == len(titles) + else: + 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): + ax = fig.add_subplot(1, n_plots, i+1) + ax.plot(x, vals[i]) + ax.set_xlabel(r'$x$') + if n_plots > 1: + ax.set_title ( titles[i] ) + plt.show() + elif gridsize_y is not None: + x = np.linspace(0, 1, gridsize_x+1) + y = np.linspace(0, 1, gridsize_y+1) + xx, yy = np.meshgrid(x, y) + vals = [[] for fun 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: + 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() + +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 _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 From 02e7cc5e92d4710f7c425bf465f68a479c7988f8 Mon Sep 17 00:00:00 2001 From: jowezarek Date: Mon, 24 Mar 2025 15:40:17 +0100 Subject: [PATCH 04/16] fix import typo --- exercises/poisson_2d_I.ipynb | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/exercises/poisson_2d_I.ipynb b/exercises/poisson_2d_I.ipynb index 0c2ec2d..04d79ce 100644 --- a/exercises/poisson_2d_I.ipynb +++ b/exercises/poisson_2d_I.ipynb @@ -130,7 +130,7 @@ "metadata": {}, "outputs": [], "source": [ - "from utils import H1BoundaryProjector2D_2\n", + "from utils import H1BoundaryProjector2D\n", "from psydac.linalg.basic import IdentityOperator\n", "\n", "P_0 = H1BoundaryProjector2D_2(V, V_h.vector_space)\n", @@ -307,5 +307,5 @@ ], "metadata": {}, "nbformat": 4, - "nbformat_minor": 2 + "nbformat_minor": 4 } From f3475347f1e9535da796c577b226371e4c94ffd9 Mon Sep 17 00:00:00 2001 From: jowezarek Date: Mon, 24 Mar 2025 15:44:30 +0100 Subject: [PATCH 05/16] fix import type 2 --- exercises/poisson_2d_I.ipynb | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/exercises/poisson_2d_I.ipynb b/exercises/poisson_2d_I.ipynb index 04d79ce..6a95fee 100644 --- a/exercises/poisson_2d_I.ipynb +++ b/exercises/poisson_2d_I.ipynb @@ -133,7 +133,7 @@ "from utils import H1BoundaryProjector2D\n", "from psydac.linalg.basic import IdentityOperator\n", "\n", - "P_0 = H1BoundaryProjector2D_2(V, V_h.vector_space)\n", + "P_0 = H1BoundaryProjector2D(V, V_h.vector_space)\n", "\n", "I_0 = IdentityOperator(V_h.vector_space)\n", "P_Gamma = I_0 - P_0" From bdfac397fd20609e63a0680ba1279bda48ecc58f Mon Sep 17 00:00:00 2001 From: jowezarek Date: Tue, 25 Mar 2025 14:28:50 +0100 Subject: [PATCH 06/16] minor beauty updates --- exercises/poisson_2d_I.ipynb | 70 ++++++++------- exercises/utils.py | 166 ++++++++++++++++------------------- 2 files changed, 114 insertions(+), 122 deletions(-) diff --git a/exercises/poisson_2d_I.ipynb b/exercises/poisson_2d_I.ipynb index 6a95fee..498df0b 100644 --- a/exercises/poisson_2d_I.ipynb +++ b/exercises/poisson_2d_I.ipynb @@ -4,7 +4,7 @@ "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.**" + "**Note: This notebook must be changed by the student. As it is now, it does not approximate the solution to the Poisson problem!**" ] }, { @@ -13,21 +13,25 @@ "source": [ "# 2D Poisson Solver on the Unit Square\n", "\n", - "In the exercise we write a solver for the 2D Poisson problem on the unit square using Psydac's bilinear form interface.\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", + "\\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", + " \\phi &= 0 \\quad \\text{ on }\\partial\\Omega,\n", + "\\end{align*}\n", + "for the specific choice\n", + "\\begin{align*}\n", + " f(x, y) = 8\\pi^2\\sin(2\\pi x)\\sin(2\\pi y)\n", + "\\end{align*}\n", "\n", "## The Variational Formulation\n", "\n", "The corresponding variational formulation reads\n", "\n", - "\\begin{align}\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", + "\\end{align*}\n", "\n", "where \n", "- $a(\\phi,\\psi) \\coloneqq \\int_{\\Omega} \\nabla\\phi\\cdot\\nabla\\psi ~ d\\Omega$,\n", @@ -52,19 +56,18 @@ "from sympde.topology import elements_of, Square\n", "from sympde.topology import ScalarFunctionSpace\n", "\n", - "from sympy import pi, sin\n", - "\n", "domain = Square('S', bounds1=(0, 1), bounds2=(0, 1))\n", "\n", - "V = ScalarFunctionSpace('V', domain)\n", - "x,y = domain.coordinates\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", + "a = BilinearForm((phi, psi), integral(domain, dot(grad(phi), grad(psi))))\n", "\n", "# linear form\n", + "from sympy import pi, sin\n", + "x,y = domain.coordinates\n", "f = 8 * (pi**2) * sin(2 * pi * x) * sin(2 * pi * y)\n", "\n", "l = LinearForm(psi, integral(domain, f*psi))" @@ -85,7 +88,10 @@ "metadata": {}, "outputs": [], "source": [ - "from psydac.api.discretization import discretize" + "from psydac.api.discretization import discretize\n", + "from psydac.api.settings import PSYDAC_BACKEND_GPYCCEL\n", + "\n", + "backend = PSYDAC_BACKEND_GPYCCEL" ] }, { @@ -94,8 +100,8 @@ "metadata": {}, "outputs": [], "source": [ - "ncells = [8, 8] # Bspline cells\n", - "degree = [2, 2] # Bspline degree" + "ncells = [16, 16] # Bspline cells\n", + "degree = [3, 3] # Bspline degree" ] }, { @@ -104,10 +110,6 @@ "metadata": {}, "outputs": [], "source": [ - "from psydac.api.settings import PSYDAC_BACKEND_GPYCCEL\n", - "\n", - "backend = PSYDAC_BACKEND_GPYCCEL\n", - "\n", "domain_h = discretize(domain, ncells=ncells, periodic=[False, False])\n", "V_h = discretize(V, domain_h, degree=degree)\n", "\n", @@ -201,10 +203,12 @@ "metadata": {}, "outputs": [], "source": [ - "from sympde.expr import find, EssentialBC, Norm, SemiNorm\n", - "from psydac.fem.basic import FemField\n", + "from psydac.fem.basic import FemField\n", + "\n", + "from sympde.expr import Norm, SemiNorm\n", "\n", "phi_ex = sin(2*pi*x) * sin(2*pi*y)\n", + "phi_h_FemField = FemField(V_h, phi_h)\n", "\n", "error = phi_ex - phi\n", "\n", @@ -212,10 +216,10 @@ "l2norm = Norm(error, domain, kind='l2')\n", "\n", "# discretize the norm\n", - "l2norm_h = discretize(l2norm, domain_h, V_h)\n", + "l2norm_h = discretize(l2norm, domain_h, V_h, backend=backend)\n", "\n", "# assemble the norm\n", - "l2_error = l2norm_h.assemble(phi=FemField(V_h, phi_h))\n", + "l2_error = l2norm_h.assemble(phi=phi_h_FemField)\n", "\n", "# print the result\n", "print(l2_error)" @@ -241,7 +245,7 @@ "h1norm_h = discretize(h1norm, domain_h, V_h)\n", "\n", "# assemble the norm\n", - "h1_error = h1norm_h.assemble(phi=FemField(V_h, phi_h))\n", + "h1_error = h1norm_h.assemble(phi=phi_h_FemField)\n", "\n", "# print the result\n", "print(h1_error)" @@ -267,7 +271,7 @@ "h1norm_h = discretize(h1norm, domain_h, V_h)\n", "\n", "# assemble the norm\n", - "h1_error = h1norm_h.assemble(phi=FemField(V_h, phi_h))\n", + "h1_error = h1norm_h.assemble(phi=phi_h_FemField)\n", "\n", "# print the result\n", "print(h1_error)" @@ -288,19 +292,19 @@ "metadata": {}, "outputs": [], "source": [ - "from utils import plot\n", "from sympy import lambdify\n", "\n", + "from utils import plot\n", + "\n", "phi_ex_fun = lambdify(domain.coordinates, phi_ex)\n", - "phi_h_fun = FemField(V_h, phi_h)\n", - "error = lambda x, y: abs(phi_ex_fun(x, y) - phi_h_fun(x, y))\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_fun, 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", + " 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", ")" ] } diff --git a/exercises/utils.py b/exercises/utils.py index c777105..5123d2a 100644 --- a/exercises/utils.py +++ b/exercises/utils.py @@ -12,10 +12,17 @@ from psydac.linalg.stencil import StencilVectorSpace -def get_unit_vector(v, n1, n2, pads1, pads2): +def get_unit_vector(v, n1, n2, n3, pads1, pads2, pads3): v *= 0.0 - v._data[pads1+n1, pads2+n2] = 1. + 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 @@ -42,97 +49,23 @@ def toarray(A): for k, Wk in enumerate(codomain_blocks): w *= 0. v = Wk.zeros() - 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, pads1, pads2) - 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() + 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 -def plot(gridsize_x, title, funs, titles, gridsize_y=None, surface_plot=False): - - if gridsize_y is None: - x = np.linspace(0, 1, gridsize_x+1) - vals = [[] for fun in funs] - for i, fun in enumerate(funs): - for xi in x: - vals[i].append(fun(xi)) - - n_plots = len(funs) - if n_plots > 1: - assert n_plots == len(titles) - else: - 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): - ax = fig.add_subplot(1, n_plots, i+1) - ax.plot(x, vals[i]) - ax.set_xlabel(r'$x$') - if n_plots > 1: - ax.set_title ( titles[i] ) - plt.show() - elif gridsize_y is not None: - x = np.linspace(0, 1, gridsize_x+1) - y = np.linspace(0, 1, gridsize_y+1) - xx, yy = np.meshgrid(x, y) - vals = [[] for fun 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: - 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() - class H1BoundaryProjector2D(LinearOperator): def __init__(self, V0, V0_vs, periodic=[False, False]): @@ -145,6 +78,8 @@ def __init__(self, V0, V0_vs, periodic=[False, False]): self._BC = self._get_BC() + def copy(self): + return H1BoundaryProjector2D(self._space, self._domain, self._periodic) def _get_BC(self): @@ -195,3 +130,56 @@ def dot(self, v, out=None): v.copy(out=out) apply_essential_bc(out, *BC) 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() From 31a2c12da583228b9956a120b38f2735ff94da71 Mon Sep 17 00:00:00 2001 From: jowezarek Date: Mon, 31 Mar 2025 14:38:29 +0200 Subject: [PATCH 07/16] add sheet1 exercise 1.4 --- exercises/time_harmonic_maxwell_2d.ipynb | 290 +++++++++++++++++++++++ exercises/utils.py | 75 +++++- 2 files changed, 363 insertions(+), 2 deletions(-) create mode 100644 exercises/time_harmonic_maxwell_2d.ipynb diff --git a/exercises/time_harmonic_maxwell_2d.ipynb b/exercises/time_harmonic_maxwell_2d.ipynb new file mode 100644 index 0000000..1b1176c --- /dev/null +++ b/exercises/time_harmonic_maxwell_2d.ipynb @@ -0,0 +1,290 @@ +{ + "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", + "\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_0 = HcurlBoundaryProjector2D(V, V_h.vector_space)\n", + "\n", + "I_0 = IdentityOperator(V_h.vector_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", + "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/utils.py b/exercises/utils.py index 5123d2a..868d833 100644 --- a/exercises/utils.py +++ b/exercises/utils.py @@ -8,8 +8,8 @@ from psydac.api.essential_bc import apply_essential_bc from psydac.linalg.basic import LinearOperator, Vector -from psydac.linalg.block import BlockVectorSpace -from psydac.linalg.stencil import StencilVectorSpace +from psydac.linalg.block import BlockVectorSpace, BlockLinearOperator +from psydac.linalg.stencil import StencilVectorSpace, StencilMatrix def get_unit_vector(v, n1, n2, n3, pads1, pads2, pads3): @@ -131,6 +131,77 @@ def dot(self, v, out=None): 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) From cada87e016a7d5cbc0cb58289ed136cb0a4d9ec2 Mon Sep 17 00:00:00 2001 From: jowezarek Date: Mon, 31 Mar 2025 14:39:24 +0200 Subject: [PATCH 08/16] modify rhs of sheet1 exercise 1.3 (prev. exercise 2.2) --- exercises/poisson_2d_I.ipynb | 41 +++++++++++++++++++++--------------- 1 file changed, 24 insertions(+), 17 deletions(-) diff --git a/exercises/poisson_2d_I.ipynb b/exercises/poisson_2d_I.ipynb index 498df0b..fa93eb8 100644 --- a/exercises/poisson_2d_I.ipynb +++ b/exercises/poisson_2d_I.ipynb @@ -21,7 +21,7 @@ "\\end{align*}\n", "for the specific choice\n", "\\begin{align*}\n", - " f(x, y) = 8\\pi^2\\sin(2\\pi x)\\sin(2\\pi y)\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", @@ -63,12 +63,12 @@ "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", + "a = BilinearForm((phi, psi), integral(domain, phi*psi))\n", "\n", "# linear form\n", - "from sympy import pi, sin\n", + "from sympy import pi, sin, cos, exp\n", "x,y = domain.coordinates\n", - "f = 8 * (pi**2) * sin(2 * pi * x) * sin(2 * pi * y)\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))" ] @@ -167,12 +167,18 @@ "metadata": {}, "outputs": [], "source": [ - "from psydac.linalg.solvers import inverse\n", + "import time\n", + "\n", + "from psydac.linalg.solvers import inverse\n", "\n", "tol = 1e-12\n", "maxiter = 1000\n", "\n", - "phi_h = inverse(A_bc, 'cg', tol=tol, maxiter=maxiter) @ f_bc" + "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()" ] }, { @@ -186,7 +192,7 @@ "In this example, the analytical solution is given by\n", "\n", "$$\n", - "phi_{ex}(x, y) = \\sin(2\\pi x)\\sin(2\\pi y)\n", + "phi_{ex}(x, y) = x(x-1)\\sin(\\pi y)e^{xy}\n", "$$" ] }, @@ -207,7 +213,7 @@ "\n", "from sympde.expr import Norm, SemiNorm\n", "\n", - "phi_ex = sin(2*pi*x) * sin(2*pi*y)\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", @@ -219,10 +225,7 @@ "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)\n", - "\n", - "# print the result\n", - "print(l2_error)" + "l2_error = l2norm_h.assemble(phi=phi_h_FemField)" ] }, { @@ -245,10 +248,7 @@ "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(h1_error)" + "h1semi_error = h1norm_h.assemble(phi=phi_h_FemField)" ] }, { @@ -274,7 +274,14 @@ "h1_error = h1norm_h.assemble(phi=phi_h_FemField)\n", "\n", "# print the result\n", - "print(h1_error)" + "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 ) )" ] }, { From 6dcdb41017a4acd439d3c2735a400ee6e41233b3 Mon Sep 17 00:00:00 2001 From: jowezarek Date: Mon, 31 Mar 2025 15:42:57 +0200 Subject: [PATCH 09/16] add sheet2 exercise 2.3 --- ...ipynb => time_harmonic_maxwell_2d_I.ipynb} | 17 +- exercises/time_harmonic_maxwell_2d_II.ipynb | 299 ++++++++++++++++++ 2 files changed, 308 insertions(+), 8 deletions(-) rename exercises/{time_harmonic_maxwell_2d.ipynb => time_harmonic_maxwell_2d_I.ipynb} (94%) create mode 100644 exercises/time_harmonic_maxwell_2d_II.ipynb diff --git a/exercises/time_harmonic_maxwell_2d.ipynb b/exercises/time_harmonic_maxwell_2d_I.ipynb similarity index 94% rename from exercises/time_harmonic_maxwell_2d.ipynb rename to exercises/time_harmonic_maxwell_2d_I.ipynb index 1b1176c..2512ba8 100644 --- a/exercises/time_harmonic_maxwell_2d.ipynb +++ b/exercises/time_harmonic_maxwell_2d_I.ipynb @@ -12,6 +12,7 @@ "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", @@ -100,7 +101,7 @@ "from psydac.api.discretization import discretize\n", "from psydac.api.settings import PSYDAC_BACKEND_GPYCCEL\n", "\n", - "backend = PSYDAC_BACKEND_GPYCCEL" + "backend = PSYDAC_BACKEND_GPYCCEL" ] }, { @@ -180,7 +181,7 @@ "\n", "from psydac.linalg.solvers import inverse\n", "\n", - "tol = 1e-12\n", + "tol = 1e-9\n", "maxiter = 1000\n", "\n", "A_bc_inv = inverse(A_bc, 'cg', tol=tol, maxiter=maxiter)\n", @@ -217,19 +218,19 @@ "\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", + "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", + "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", + "l2norm = Norm(error, domain, kind='l2')\n", "\n", "# discretize the norm\n", - "l2norm_h = discretize(l2norm, domain_h, V_h, backend=backend)\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", + "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", diff --git a/exercises/time_harmonic_maxwell_2d_II.ipynb b/exercises/time_harmonic_maxwell_2d_II.ipynb new file mode 100644 index 0000000..caae8bd --- /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_0 = HcurlBoundaryProjector2D(V1, V1_h.vector_space)\n", + "\n", + "I_0 = IdentityOperator(V1_h.vector_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-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 +} From cb5f6343fbba1a482ec75cd5554d1c98f5189da7 Mon Sep 17 00:00:00 2001 From: jowezarek Date: Mon, 31 Mar 2025 16:12:31 +0200 Subject: [PATCH 10/16] add sheet2 exercise 2.2 --- exercises/poisson_2d_I.ipynb | 4 +- exercises/poisson_2d_II.ipynb | 329 ++++++++++++++++++++ exercises/time_harmonic_maxwell_2d_I.ipynb | 4 +- exercises/time_harmonic_maxwell_2d_II.ipynb | 4 +- 4 files changed, 335 insertions(+), 6 deletions(-) create mode 100644 exercises/poisson_2d_II.ipynb diff --git a/exercises/poisson_2d_I.ipynb b/exercises/poisson_2d_I.ipynb index fa93eb8..c6f704e 100644 --- a/exercises/poisson_2d_I.ipynb +++ b/exercises/poisson_2d_I.ipynb @@ -12,6 +12,7 @@ "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", @@ -53,8 +54,7 @@ "source": [ "from sympde.calculus import dot, grad\n", "from sympde.expr import BilinearForm, LinearForm, integral\n", - "from sympde.topology import elements_of, Square\n", - "from sympde.topology import ScalarFunctionSpace\n", + "from sympde.topology import elements_of, Square, ScalarFunctionSpace\n", "\n", "domain = Square('S', bounds1=(0, 1), bounds2=(0, 1))\n", "\n", diff --git a/exercises/poisson_2d_II.ipynb b/exercises/poisson_2d_II.ipynb new file mode 100644 index 0000000..444fb88 --- /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.vector_space)\n", + "\n", + "I_0 = IdentityOperator(V0_h.vector_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/time_harmonic_maxwell_2d_I.ipynb b/exercises/time_harmonic_maxwell_2d_I.ipynb index 2512ba8..eb228e2 100644 --- a/exercises/time_harmonic_maxwell_2d_I.ipynb +++ b/exercises/time_harmonic_maxwell_2d_I.ipynb @@ -61,9 +61,9 @@ "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", + "domain = Square('S', bounds1=(0, 1), bounds2=(0, 1))\n", "\n", - "V = VectorFunctionSpace('V', domain, kind='hcurl')\n", + "V = VectorFunctionSpace('V', domain, kind='hcurl')\n", "\n", "E, v = elements_of(V, names='E, v')\n", "\n", diff --git a/exercises/time_harmonic_maxwell_2d_II.ipynb b/exercises/time_harmonic_maxwell_2d_II.ipynb index caae8bd..ec7a565 100644 --- a/exercises/time_harmonic_maxwell_2d_II.ipynb +++ b/exercises/time_harmonic_maxwell_2d_II.ipynb @@ -61,8 +61,8 @@ "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", + "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", From e92e5ec8b1aae501c17710e54ad634115f73ffeb Mon Sep 17 00:00:00 2001 From: jowezarek Date: Tue, 1 Apr 2025 11:10:40 +0200 Subject: [PATCH 11/16] add sheet3 exercise 3.2 --- exercises/electro_static_2d.ipynb | 281 ++++++++++++++++++++ exercises/time_harmonic_maxwell_2d_I.ipynb | 10 +- exercises/time_harmonic_maxwell_2d_II.ipynb | 10 +- exercises/utils.py | 95 ++++++- 4 files changed, 382 insertions(+), 14 deletions(-) create mode 100644 exercises/electro_static_2d.ipynb diff --git a/exercises/electro_static_2d.ipynb b/exercises/electro_static_2d.ipynb new file mode 100644 index 0000000..15afaa7 --- /dev/null +++ b/exercises/electro_static_2d.ipynb @@ -0,0 +1,281 @@ +{ + "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.vector_space)\n", + "#P_Hcurl = HcurlBoundaryProjector2D(V1, V1_h.vector_space)\n", + "\n", + "#I0 = IdentityOperator(V0_h.vector_space)\n", + "#I1 = IdentityOperator(V1_h.vector_space)\n", + "\n", + "#P_H1_Gamma = I0 - P_H1\n", + "#P_Hcurl_Gamma = I1 - P_Hcurl\n", + "\n", + "# Modified Operators for Projection Method (using P_H1, P_Hcurl, P_H1_Gamma, P_Hcurl_Gamma)\n", + "# <-\n", + "# ...\n", + "# <-\n", + "\n", + "# Inverse M0 mass matrix\n", + "from psydac.linalg.solvers import inverse\n", + "# <-\n", + "\n", + "# System Matrix A\n", + "#A = # <-\n", + "\n", + "# rhs vector f\n", + "#rho_coeffs = P0(rho_lambdified).coeffs\n", + "#f = # <-" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import time\n", + "\n", + "tol = 1e-10\n", + "maxiter = 1000\n", + "\n", + "from utils import get_M1_block_kron_solver_2D\n", + "#pc = get_M1_block_kron_solver_2D(V1_h.vector_space, ncells, degree, periodic=[False, False])\n", + "#M1_0_inv = inverse(M1_0, 'pcg', pc=pc, tol=1e-11, maxiter=1000)\n", + "#A_inv = inverse(A, 'pcg', pc=M1_0_inv, tol=tol, maxiter=maxiter)\n", + "\n", + "t0 = time.time()\n", + "#E_h = A_inv @ f\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_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", + "if False:\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_I.ipynb b/exercises/time_harmonic_maxwell_2d_I.ipynb index eb228e2..5d2c5a6 100644 --- a/exercises/time_harmonic_maxwell_2d_I.ipynb +++ b/exercises/time_harmonic_maxwell_2d_I.ipynb @@ -145,10 +145,10 @@ "from utils import HcurlBoundaryProjector2D\n", "from psydac.linalg.basic import IdentityOperator\n", "\n", - "P_0 = HcurlBoundaryProjector2D(V, V_h.vector_space)\n", + "P_Hcurl = HcurlBoundaryProjector2D(V, V_h.vector_space)\n", "\n", - "I_0 = IdentityOperator(V_h.vector_space)\n", - "P_Gamma = I_0 - P_0" + "I1 = IdentityOperator(V_h.vector_space)\n", + "P_Gamma = I1 - P_Hcurl" ] }, { @@ -160,8 +160,8 @@ "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" + "A_bc = P_Hcurl @ A @ P_Hcurl + P_Gamma\n", + "f_bc = P_Hcurl @ f" ] }, { diff --git a/exercises/time_harmonic_maxwell_2d_II.ipynb b/exercises/time_harmonic_maxwell_2d_II.ipynb index ec7a565..2a84277 100644 --- a/exercises/time_harmonic_maxwell_2d_II.ipynb +++ b/exercises/time_harmonic_maxwell_2d_II.ipynb @@ -156,10 +156,10 @@ "from utils import HcurlBoundaryProjector2D\n", "from psydac.linalg.basic import IdentityOperator\n", "\n", - "P_0 = HcurlBoundaryProjector2D(V1, V1_h.vector_space)\n", + "P_Hcurl = HcurlBoundaryProjector2D(V1, V1_h.vector_space)\n", "\n", - "I_0 = IdentityOperator(V1_h.vector_space)\n", - "P_Gamma = I_0 - P_0" + "I1 = IdentityOperator(V1_h.vector_space)\n", + "P_Gamma = I1 - P_Hcurl" ] }, { @@ -168,8 +168,8 @@ "metadata": {}, "outputs": [], "source": [ - "A_bc = P_0 @ A @ P_0 + P_Gamma\n", - "f_bc = P_0 @ f" + "A_bc = P_Hcurl @ A @ P_Hcurl + P_Gamma\n", + "f_bc = P_Hcurl @ f" ] }, { diff --git a/exercises/utils.py b/exercises/utils.py index 868d833..27c4068 100644 --- a/exercises/utils.py +++ b/exercises/utils.py @@ -1,15 +1,20 @@ -import h5py import numpy as np import matplotlib.pyplot as plt from matplotlib import cm, colors -from sympde.expr import EssentialBC -from sympde.topology import element_of +from scipy.sparse import dia_matrix +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.stencil import StencilVectorSpace, StencilMatrix +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): @@ -254,3 +259,85 @@ def plot(gridsize_x, gridsize_y, title, funs, titles, surface_plot=False): 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 From f2f6b200afeb80ea227b01e495d0e72d20640408 Mon Sep 17 00:00:00 2001 From: jowezarek Date: Wed, 2 Apr 2025 16:30:35 +0200 Subject: [PATCH 12/16] add (first version of) sheet4 exercise 4.7 --- exercises/hodge_laplace_annulus.ipynb | 432 ++++++++++++++++++++++++++ exercises/utils.py | 39 +++ 2 files changed, 471 insertions(+) create mode 100644 exercises/hodge_laplace_annulus.ipynb diff --git a/exercises/hodge_laplace_annulus.ipynb b/exercises/hodge_laplace_annulus.ipynb new file mode 100644 index 0000000..d847ada --- /dev/null +++ b/exercises/hodge_laplace_annulus.ipynb @@ -0,0 +1,432 @@ +{ + "cells": [ + { + "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\n", + "from sympde.expr import LinearForm, BilinearForm, integral\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-9)" + ] + }, + { + "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 = Ct_sp @ M2_sp @ C_sp + M1_sp @ G_sp_0 @ inv_M0_sp @ Gt_sp_0 @ M1_sp\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([[M0_sp, Gt_sp @ M1_sp, None],\n", + " [M1_sp @ G_sp, Ct_sp @ M2_sp @ C_sp, M1_sp @ hf_sp],\n", + " [None, hft_sp @ M1_sp, None]])\n", + "\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": "markdown", + "metadata": {}, + "source": [ + "# Method of manufactured solution\n", + "\n", + "In order to test our code, we employ the method of manufactured solution.\n", + "We start from the true $\\boldsymbol{A}\\in H_0(curl;\\Omega),\\ div\\ \\boldsymbol{A} = 0$\n", + "$$\n", + "\\boldsymbol{A}(x, y) = \\frac{1}{(x^2 + y^2)^{3/2}}\\left(\\begin{matrix} xy \\\\ y^2 \\end{matrix}\\right)\n", + "$$\n", + "Because the r.h.s. of the discrete variational formulation reads\n", + "$$\n", + "r.h.s. = \\left(\\begin{matrix} \\mathbb{0} \\\\ (\\boldsymbol{v},\\ \\boldsymbol{J}) \\\\ \\mathbb{0} \\end{matrix}\\right)\n", + "$$\n", + "and because $\\boldsymbol{J} = \\boldsymbol{curl}\\ curl\\ \\boldsymbol{A}$, the r.h.s. of the linear system becomes\n", + "$$\n", + "r.h.s. = \\left(\\begin{matrix} \\mathbb{0} \\\\ \\mathbb{C}^T\\mathbb{M}_2\\mathbb{C}\\mathbb{A} \\\\ \\mathbb{0} \\end{matrix}\\right)\n", + "$$\n", + "with $\\mathbb{A}$ being the coefficient vector of discrete $\\boldsymbol{A}$. If our code is correct, we will obtain approximately $\\mathbb{A}$ as second component of the solution vector." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Define exact solution A_ex for method of manufactured solution\n", + "r = lambda x, y : np.sqrt(x**2 + y**2)\n", + "A_ex_1 = lambda x, y : (x*y) / (r(x,y)**3)\n", + "A_ex_2 = lambda x, y : (y**2) / (r(x,y)**3)\n", + "\n", + "A_ex = (A_ex_1, A_ex_2)\n", + "A_ex_FemField = P1(A_ex)\n", + "A_ex_coeffs = A_ex_FemField.coeffs\n", + "\n", + "# First component of rhs\n", + "rhs1_nparray = np.zeros(V0h.nbasis)\n", + "\n", + "# Second component of rhs\n", + "rhs2_nparray = (P_Hcurl @ Ct @ M2 @ C @ A_ex_coeffs ).toarray()\n", + "\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": [ + "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", + "Ah_coeffs = array_to_psydac(Ah_nparray, V1h.coeff_space)\n", + "\n", + "error = A_ex_coeffs - Ah_coeffs\n", + "l2_error = np.sqrt( error.dot(M1 @ error) )\n", + "\n", + "# print the result\n", + "print( '> Grid :: [{ne1},{ne2}]'.format( ne1=ncells[0], ne2=ncells[1]) )\n", + "print( '> Degree :: [{p1},{p2}]' .format( p1=degree[0], p2=degree[1] ) )\n", + "print( '> L2 error :: {:.2e}'.format( l2_error ) )\n", + "print( '' )\n", + "print( '> Solution time :: {:.3g}'.format( t1-t0 ) )" + ] + } + ], + "metadata": {}, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/exercises/utils.py b/exercises/utils.py index 27c4068..27abb27 100644 --- a/exercises/utils.py +++ b/exercises/utils.py @@ -3,6 +3,7 @@ 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 @@ -341,3 +342,41 @@ def get_M1_block_kron_solver_2D(V1, ncells, degree, periodic): (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 From f9537c01b7f07ca100a803f2429a28ec06b3fe59 Mon Sep 17 00:00:00 2001 From: jowezarek Date: Thu, 3 Apr 2025 16:48:18 +0200 Subject: [PATCH 13/16] complete sheet4 exercise 4.7 --- exercises/hodge_laplace_annulus.ipynb | 229 +++++++++++++++++++------- 1 file changed, 168 insertions(+), 61 deletions(-) diff --git a/exercises/hodge_laplace_annulus.ipynb b/exercises/hodge_laplace_annulus.ipynb index d847ada..2009457 100644 --- a/exercises/hodge_laplace_annulus.ipynb +++ b/exercises/hodge_laplace_annulus.ipynb @@ -1,5 +1,71 @@ { "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "**Note: This notebook must be changed by the student. As it is now, it does not approximate the solution to the Poisson problem!**" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# 2D Solver for the Magneto-Static Problem on an Annulus\n", + "\n", + "In this exercise we write a solver for the 2D magneto-static problem on an annulus (r_min=.3, r_max=1.) $\\Omega$.\n", + "\n", + "\\begin{align*}\n", + " &\\text{Given }\\boldsymbol{J}\\in H_0(curl;\\Omega)\\cap(\\mathfrak{H}^1)^{\\perp}\\text{ s.th. }div\\ \\boldsymbol{J}=0, \\\\\n", + " &\\text{find }B\\in H(\\boldsymbol{curl};\\Omega)\\text{ satisfying} \\\\\n", + " &\\boldsymbol{curl}\\ B = \\boldsymbol{J} \\quad \\text{ in }\\Omega,\n", + "\\end{align*}\n", + "\n", + "for the specific choice\n", + "\n", + "\\begin{align*}\n", + " r(x, y) &= \\sqrt{x^2 + y^2}, \\\\\n", + " B(x, y) &= \\frac{-x}{r^3}, \\\\\n", + " \\boldsymbol{J}(x, y) &= \\boldsymbol{curl}\\ B.\n", + "\\end{align*}\n", + "\n", + "Note that $B$ is chosen to be the $curl$ of $\\boldsymbol{A}$ such that $\\boldsymbol{A}\\in H_0(curl;\\Omega)$ and $div\\ \\boldsymbol{A} = 0$ for\n", + "\\begin{align}\n", + " \\boldsymbol{A}(x, y) = \\frac{1}{r^3}\\left(\\begin{matrix} xy \\\\ y^2 \\end{matrix}\\right).\n", + "\\end{align} \n", + "\n", + "This problem is ill-posed. We have shown, that the following auxiliary problem is well-posed, and that its solution $\\boldsymbol{A}$ satisfies $\\boldsymbol{curl}\\ curl\\ \\boldsymbol{A} = \\boldsymbol{J}$, hence $B = curl\\ \\boldsymbol{A}$ is a solution to the original problem.\n", + "\n", + "\\begin{align*}\n", + " &\\text{Given }\\boldsymbol{J}\\in H_0(curl;\\Omega)\\cap(\\mathfrak{H}^1)^{\\perp}\\text{ s.th. }div\\ \\boldsymbol{J}=0, \\\\\n", + " &\\text{find }\\boldsymbol{A}\\in D(\\mathcal{L}^1)\\cap(\\mathfrak{H}^1)^{\\perp}\\text{ satisfying} \\\\\n", + " &\\mathcal{L}^1(\\boldsymbol{A}) = \\boldsymbol{J} \\quad \\text{ in }\\Omega.\n", + "\\end{align*}\n", + "\n", + "## The Discrete Variational Formulation\n", + "\n", + "The corresponding discrete variational formulation reads\n", + "\n", + "\\begin{align*}\n", + " &\\quad\\text{For }\\boldsymbol{J}\\in V_h^1\\cap(\\mathfrak{H}_h^1)^{\\perp}\\text{ s.th. }\\widetilde{div}_h\\ \\boldsymbol{J}=0, \\\\\n", + " &\\quad\\text{find } (\\sigma,\\boldsymbol{A},\\boldsymbol{p})\\in V_h^0\\times V_h^1\\times\\mathfrak{H}_h^1 \\text{ satisfying} \\\\\n", + " (\\textbf{HL1}_h)\n", + " &\\begin{cases}\n", + " \\quad\\quad\\quad\\ \\ \\langle\\sigma,\\tau\\rangle - \\quad\\ \\ \\langle\\boldsymbol{A},\\boldsymbol{grad}\\ \\tau\\rangle &= 0 \\ \\quad\\quad\\forall\\tau\\in V_h^0, \\\\\n", + " -\\langle\\boldsymbol{grad}\\ \\sigma,\\boldsymbol{v}\\rangle + \\langle curl\\ \\boldsymbol{A},curl\\ \\boldsymbol{v}\\rangle + \\langle\\boldsymbol{p},\\boldsymbol{v}\\rangle &= \\langle\\boldsymbol{J},\\boldsymbol{v}\\rangle \\ \\forall\\boldsymbol{v}\\in V_h^1, \\\\\n", + " \\quad\\quad\\quad\\quad\\quad\\quad\\quad\\quad\\quad\\quad\\quad\\langle\\boldsymbol{A},\\boldsymbol{q}\\rangle &= 0 \\ \\quad\\quad\\forall \\boldsymbol{q}\\in\\mathfrak{H}_h^1.\n", + " \\end{cases}\n", + "\\end{align*}\n", + "\n", + "The discrete solution $\\boldsymbol{A}$ will then satisfy $\\widetilde{\\boldsymbol{curl}}_h\\ curl\\ \\boldsymbol{A} = \\boldsymbol{J}$, hence $B = curl\\ \\boldsymbol{A}$ will satisfy $\\widetilde{\\boldsymbol{curl}}_h B = \\boldsymbol{J}$.\n", + "\n", + "## Method of manufactured solution\n", + "\n", + "To verify that our code works, we employ the method of manufactured solutions. Given the vector field $\\boldsymbol{A}_{ex}$, see (1), we compute the r.h.s. via $\\langle\\boldsymbol{J},\\ \\boldsymbol{v}\\rangle = \\langle\\widetilde{\\boldsymbol{curl}}_h\\ curl\\ \\boldsymbol{A}_{ex},\\ \\boldsymbol{v}\\rangle = \\langle curl\\ \\boldsymbol{A}_{ex},\\ curl\\ \\boldsymbol{v}\\rangle$.\n", + "\n", + "We then report the $L^2$ errors $\\| \\boldsymbol{A}_h - \\boldsymbol{A}_{ex} \\|$ and $\\| B_h - B_{ex} \\|$." + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -46,8 +112,8 @@ "metadata": {}, "outputs": [], "source": [ - "from sympde.calculus import dot\n", - "from sympde.expr import LinearForm, BilinearForm, integral\n", + "from sympde.calculus import dot, curl\n", + "from sympde.expr import LinearForm, BilinearForm, integral, Norm\n", "from sympde.topology import elements_of, Derham\n", "\n", "derham = Derham(domain, sequence=['h1', 'hcurl', 'l2'])\n", @@ -168,7 +234,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Boundary Conditions and inverse Matrices\n", + "## Boundary Conditions and Inverse Matrices\n", "\n", "Let us denote the interior coefficients by a subscript $I$ and the exterior coefficients by a subscript $\\Gamma$. Suppose we want to solve the linear system\n", "$$\n", @@ -208,7 +274,7 @@ "from psydac.linalg.solvers import inverse\n", "\n", "M0_0 = P_H1 @ M0 @ P_H1 + P_H1_Gamma\n", - "M0_inv = inverse(M0_0, 'cg', maxiter=1000, tol=1e-9)" + "M0_inv = inverse(M0_0, 'cg', maxiter=1000, tol=1e-12)" ] }, { @@ -220,8 +286,8 @@ "Psydac does not offer a direct way to access eigenvalues and -vectors of matrices.\n", "However, we can transform our operators to scipy.sparse format and use scipy's `eigsh` method.\n", "\n", - "We want our solution vector field $A$ to be an element of $(\\mathfrak{h}_h^1)^{\\perp}$, \n", - "and it holds $\\ker\\mathcal{L}_h^1 = \\mathfrak{h}_h^1$.\n", + "We want our solution vector field $A$ to be an element of $(\\mathfrak{H}_h^1)^{\\perp}$, \n", + "and it holds $\\ker\\mathcal{L}_h^1 = \\mathfrak{H}_h^1$.\n", "\n", "Because the annulus has one whole, and hence the space of harmonic forms is of dimension 1, we want to compute the first eigenvector (corresponding to the eigenvalue $0$) of the operator\n", "$$\n", @@ -271,9 +337,9 @@ "inv_M0_sp.eliminate_zeros()\n", "\n", "# For this complex operator, the projection method consists of using the modified gradient\n", - "L1_sp = Ct_sp @ M2_sp @ C_sp + M1_sp @ G_sp_0 @ inv_M0_sp @ Gt_sp_0 @ M1_sp\n", + "#L1_sp = # <- todo \n", "# ... and the projection matrices as usual\n", - "L1_sp_bc = P_Hcurl_sp @ L1_sp @ P_Hcurl_sp + P_Hcurl_Gamma_sp" + "#L1_sp_bc = P_Hcurl_sp @ L1_sp @ P_Hcurl_sp + P_Hcurl_Gamma_sp" ] }, { @@ -292,14 +358,14 @@ "from utils import get_eigenvalues\n", "\n", "dim_harmonic_space = 1\n", - "eigenvalues, eigenvectors = get_eigenvalues(dim_harmonic_space, 1e-6, L1_sp_bc, M1_sp)\n", + "#eigenvalues, eigenvectors = get_eigenvalues(dim_harmonic_space, 1e-6, L1_sp_bc, M1_sp)\n", "\n", "# A basis of the vector space of harmonic forms (hf)\n", - "hf = eigenvectors[:,0]\n", + "#hf = eigenvectors[:,0]\n", "# Obtain the eigenvector as a \"row csc_matrix (1 x n)\"\n", - "hft_sp = csc_matrix(hf)\n", + "#hft_sp = csc_matrix(hf)\n", "# ... and as a \"column csc_matrix (n x 1)\"\n", - "hf_sp = hft_sp.transpose()" + "#hf_sp = hft_sp.transpose()" ] }, { @@ -317,9 +383,9 @@ "outputs": [], "source": [ "# System matrix without taking BCs into account\n", - "A_sp = bmat([[M0_sp, Gt_sp @ M1_sp, None],\n", - " [M1_sp @ G_sp, Ct_sp @ M2_sp @ C_sp, M1_sp @ hf_sp],\n", - " [None, hft_sp @ M1_sp, None]])\n", + "#A_sp = bmat([[x, x, None], # <- todo\n", + "# [x, x, M1_sp @ x], # <- todo\n", + "# [None, x, None]]) # <- todo\n", "\n", "ZERO = csc_matrix(np.array([0]))\n", "ONE = csc_matrix(np.array([1]))\n", @@ -333,29 +399,7 @@ " [None, None, ZERO]])\n", "\n", "# System matrix taking BCs into account\n", - "A_sp_bc = P @ A_sp @ P + P_Gamma" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Method of manufactured solution\n", - "\n", - "In order to test our code, we employ the method of manufactured solution.\n", - "We start from the true $\\boldsymbol{A}\\in H_0(curl;\\Omega),\\ div\\ \\boldsymbol{A} = 0$\n", - "$$\n", - "\\boldsymbol{A}(x, y) = \\frac{1}{(x^2 + y^2)^{3/2}}\\left(\\begin{matrix} xy \\\\ y^2 \\end{matrix}\\right)\n", - "$$\n", - "Because the r.h.s. of the discrete variational formulation reads\n", - "$$\n", - "r.h.s. = \\left(\\begin{matrix} \\mathbb{0} \\\\ (\\boldsymbol{v},\\ \\boldsymbol{J}) \\\\ \\mathbb{0} \\end{matrix}\\right)\n", - "$$\n", - "and because $\\boldsymbol{J} = \\boldsymbol{curl}\\ curl\\ \\boldsymbol{A}$, the r.h.s. of the linear system becomes\n", - "$$\n", - "r.h.s. = \\left(\\begin{matrix} \\mathbb{0} \\\\ \\mathbb{C}^T\\mathbb{M}_2\\mathbb{C}\\mathbb{A} \\\\ \\mathbb{0} \\end{matrix}\\right)\n", - "$$\n", - "with $\\mathbb{A}$ being the coefficient vector of discrete $\\boldsymbol{A}$. If our code is correct, we will obtain approximately $\\mathbb{A}$ as second component of the solution vector." + "#A_sp_bc = P @ A_sp @ P + P_Gamma" ] }, { @@ -364,32 +408,33 @@ "metadata": {}, "outputs": [], "source": [ - "# Define exact solution A_ex for method of manufactured solution\n", - "r = lambda x, y : np.sqrt(x**2 + y**2)\n", - "A_ex_1 = lambda x, y : (x*y) / (r(x,y)**3)\n", - "A_ex_2 = lambda x, y : (y**2) / (r(x,y)**3)\n", + "from sympy import sqrt, Tuple, Matrix\n", + "x, y = domain.coordinates\n", + "r_sympy = sqrt(x**2 + y**2)\n", + "A_ex_1_sympy = (x*y) / (r_sympy**3)\n", + "A_ex_2_sympy = (y**2) / (r_sympy**3)\n", + "A_ex_sympy = Tuple(A_ex_1_sympy, A_ex_2_sympy)\n", "\n", - "A_ex = (A_ex_1, A_ex_2)\n", - "A_ex_FemField = P1(A_ex)\n", - "A_ex_coeffs = A_ex_FemField.coeffs\n", + "#l = # <- todo\n", + "#lh = discretize(l, domain_h, V1h)\n", "\n", "# First component of rhs\n", "rhs1_nparray = np.zeros(V0h.nbasis)\n", "\n", "# Second component of rhs\n", - "rhs2_nparray = (P_Hcurl @ Ct @ M2 @ C @ A_ex_coeffs ).toarray()\n", + "#rhs2_nparray = (P_Hcurl @ lh.assemble()).toarray()\n", "\n", "# Third component of rhs\n", "rhs3_nparray = np.zeros(dim_harmonic_space)\n", "\n", "# Full rhs\n", - "rhs_nparray = np.block([rhs1_nparray, rhs2_nparray, rhs3_nparray])\n", + "#rhs_nparray = np.block([rhs1_nparray, rhs2_nparray, rhs3_nparray])\n", "\n", "# -------------- direct solve with scipy spsolve ---------------\n", "import time\n", "from scipy.sparse.linalg import spsolve\n", "t0 = time.time()\n", - "sol_nparray = spsolve(A_sp_bc.asformat('csr'), rhs_nparray)\n", + "#sol_nparray = spsolve(A_sp_bc.asformat('csr'), rhs_nparray)\n", "t1 = time.time()\n", "# --------------------------------------------------------------" ] @@ -407,22 +452,84 @@ "metadata": {}, "outputs": [], "source": [ - "sigmah_nparray = sol_nparray[:V0h.nbasis]\n", - "Ah_nparray = sol_nparray[V0h.nbasis:V0h.nbasis + V1h.nbasis]\n", - "ph_nparray = sol_nparray[V0h.nbasis+V1h.nbasis:]\n", + "# print the result\n", + "print( '> Grid :: [{ne1},{ne2}]'.format( ne1=ncells[0], ne2=ncells[1]) )\n", + "print( '> Degree :: [{p1},{p2}]' .format( p1=degree[0], p2=degree[1] ) )\n", + "print( '> Solution time :: {:.3g}'.format( t1-t0 ) )\n", + "print( '' )\n", "\n", - "from psydac.linalg.utilities import array_to_psydac\n", - "Ah_coeffs = array_to_psydac(Ah_nparray, V1h.coeff_space)\n", + "#sigmah_nparray = sol_nparray[:V0h.nbasis]\n", + "#Ah_nparray = sol_nparray[V0h.nbasis:V0h.nbasis + V1h.nbasis]\n", + "#ph_nparray = sol_nparray[V0h.nbasis+V1h.nbasis:]\n", + "\n", + "from psydac.linalg.utilities import array_to_psydac\n", + "from psydac.fem.basic import FemField\n", + "#Ah_coeffs = array_to_psydac(Ah_nparray, V1h.coeff_space)\n", + "\n", + "# ----- Ah vs. Aex -----\n", + "#Ah_FemField = FemField(V1h, Ah_coeffs)\n", + "error = Matrix([u1[0] - A_ex_sympy[0], u1[1] - A_ex_sympy[1]])\n", + "l2norm = Norm(error, domain, kind='l2')\n", + "l2norm_h = discretize(l2norm, domain_h, V1h, backend=backend)\n", + "#l2_error = l2norm_h.assemble(u1=Ah_FemField)\n", + "\n", + "#print( '> || Ah - Aex || :: {:.2e}'.format( l2_error ) )\n", + "# ----------------------\n", + "\n", + "# ----- Bh vs. Bex -----\n", + "B_ex_sympy = ( -x ) / ( r_sympy**3 )\n", + "#Bh_coeffs = C @ Ah_coeffs\n", + "#Bh_FemField = FemField(V2h, Bh_coeffs)\n", + "error = u2 - B_ex_sympy\n", + "l2norm = Norm(error, domain, kind='l2')\n", + "l2norm_h = discretize(l2norm, domain_h, V2h, backend=backend)\n", + "#l2_error = l2norm_h.assemble(u2=Bh_FemField)\n", + "\n", + "#print( '> || Bh - Bex || :: {:.2e}'.format( l2_error ) )\n", + "# ----------------------" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from sympy import lambdify\n", "\n", - "error = A_ex_coeffs - Ah_coeffs\n", - "l2_error = np.sqrt( error.dot(M1 @ error) )\n", + "from psydac.fem.plotting_utilities import plot_field_2d\n", "\n", - "# print the result\n", - "print( '> Grid :: [{ne1},{ne2}]'.format( ne1=ncells[0], ne2=ncells[1]) )\n", - "print( '> Degree :: [{p1},{p2}]' .format( p1=degree[0], p2=degree[1] ) )\n", - "print( '> L2 error :: {:.2e}'.format( l2_error ) )\n", - "print( '' )\n", - "print( '> Solution time :: {:.3g}'.format( t1-t0 ) )" + "Aex_lambdified = lambdify((x, y), A_ex_sympy)\n", + "#A_error_fun = lambda x, y : Aex_lambdified(x, y) - Ah_FemField(x, y)\n", + "\n", + "Bex_lambdified = lambdify((x, y), B_ex_sympy)\n", + "#B_error_fun = lambda x, y : Bex_lambdified(x, y) - Bh_FemField(x, y)\n", + "\n", + "#plot_field_2d(Ah_FemField, domain=domain, space_kind='hcurl', N_vis=250, title=r'$\\| \\boldsymbol{A}_h \\|$', filename=None, plot_type='vector_field')\n", + "#plot_field_2d(Ah_FemField, domain=domain, space_kind='hcurl', N_vis=50, title=r'$\\| \\boldsymbol{A}_h \\|$', filename=None)\n", + "\n", + "#plot_field_2d(Bh_FemField, domain=domain, space_kind='l2', N_vis=500, title=r'$\\| B_h \\|$', filename=None, surface_plot=True)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Harmonic Field\n", + "\n", + "We can also plot the harmonic field" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#hf_coeffs = array_to_psydac(hf, V1h.coeff_space)\n", + "\n", + "#plot_field_2d(stencil_coeffs=hf_coeffs, Vh=V1h, domain=domain, space_kind='hcurl', N_vis=50, title=r'$\\| \\boldsymbol{\\mathfrak{H}} \\|$', filename=None)\n", + "#plot_field_2d(stencil_coeffs=hf_coeffs, Vh=V1h, domain=domain, space_kind='hcurl', N_vis=250, title=r'$\\| \\boldsymbol{\\mathfrak{H}} \\|$', filename=None, plot_type='vector_field')" ] } ], From 763405b752fcf93a4df2a643840e83a5063b7a1f Mon Sep 17 00:00:00 2001 From: jowezarek Date: Thu, 3 Apr 2025 17:20:39 +0200 Subject: [PATCH 14/16] add sheet1 exercise 1.3 iii) --- exercises/poisson_2d_I_naturalBC.ipynb | 304 +++++++++++++++++++++++++ 1 file changed, 304 insertions(+) create mode 100644 exercises/poisson_2d_I_naturalBC.ipynb diff --git a/exercises/poisson_2d_I_naturalBC.ipynb b/exercises/poisson_2d_I_naturalBC.ipynb new file mode 100644 index 0000000..da5fabe --- /dev/null +++ b/exercises/poisson_2d_I_naturalBC.ipynb @@ -0,0 +1,304 @@ +{ + "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 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 +} From b9b53996c802618cf5a43e5434ee084509bcc9e3 Mon Sep 17 00:00:00 2001 From: jowezarek Date: Fri, 4 Apr 2025 09:52:55 +0200 Subject: [PATCH 15/16] change .vector_space to .coeff_space --- exercises/electro_static_2d.ipynb | 10 +++++----- exercises/poisson_2d_I.ipynb | 4 ++-- exercises/poisson_2d_II.ipynb | 4 ++-- exercises/poisson_2d_I_naturalBC.ipynb | 7 ------- exercises/time_harmonic_maxwell_2d_I.ipynb | 4 ++-- exercises/time_harmonic_maxwell_2d_II.ipynb | 4 ++-- 6 files changed, 13 insertions(+), 20 deletions(-) diff --git a/exercises/electro_static_2d.ipynb b/exercises/electro_static_2d.ipynb index 15afaa7..399edce 100644 --- a/exercises/electro_static_2d.ipynb +++ b/exercises/electro_static_2d.ipynb @@ -134,11 +134,11 @@ "from utils import H1BoundaryProjector2D, HcurlBoundaryProjector2D\n", "from psydac.linalg.basic import IdentityOperator\n", "\n", - "#P_H1 = H1BoundaryProjector2D(V0, V0_h.vector_space)\n", - "#P_Hcurl = HcurlBoundaryProjector2D(V1, V1_h.vector_space)\n", + "#P_H1 = H1BoundaryProjector2D(V0, V0_h.coeff_space)\n", + "#P_Hcurl = HcurlBoundaryProjector2D(V1, V1_h.coeff_space)\n", "\n", - "#I0 = IdentityOperator(V0_h.vector_space)\n", - "#I1 = IdentityOperator(V1_h.vector_space)\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", @@ -172,7 +172,7 @@ "maxiter = 1000\n", "\n", "from utils import get_M1_block_kron_solver_2D\n", - "#pc = get_M1_block_kron_solver_2D(V1_h.vector_space, ncells, degree, periodic=[False, False])\n", + "#pc = get_M1_block_kron_solver_2D(V1_h.coeff_space, ncells, degree, periodic=[False, False])\n", "#M1_0_inv = inverse(M1_0, 'pcg', pc=pc, tol=1e-11, maxiter=1000)\n", "#A_inv = inverse(A, 'pcg', pc=M1_0_inv, tol=tol, maxiter=maxiter)\n", "\n", diff --git a/exercises/poisson_2d_I.ipynb b/exercises/poisson_2d_I.ipynb index c6f704e..ca64bab 100644 --- a/exercises/poisson_2d_I.ipynb +++ b/exercises/poisson_2d_I.ipynb @@ -135,9 +135,9 @@ "from utils import H1BoundaryProjector2D\n", "from psydac.linalg.basic import IdentityOperator\n", "\n", - "P_0 = H1BoundaryProjector2D(V, V_h.vector_space)\n", + "P_0 = H1BoundaryProjector2D(V, V_h.coeff_space)\n", "\n", - "I_0 = IdentityOperator(V_h.vector_space)\n", + "I_0 = IdentityOperator(V_h.coeff_space)\n", "P_Gamma = I_0 - P_0" ] }, diff --git a/exercises/poisson_2d_II.ipynb b/exercises/poisson_2d_II.ipynb index 444fb88..e22c2d8 100644 --- a/exercises/poisson_2d_II.ipynb +++ b/exercises/poisson_2d_II.ipynb @@ -145,9 +145,9 @@ "from utils import H1BoundaryProjector2D\n", "from psydac.linalg.basic import IdentityOperator\n", "\n", - "P_0 = H1BoundaryProjector2D(V0, V0_h.vector_space)\n", + "P_0 = H1BoundaryProjector2D(V0, V0_h.coeff_space)\n", "\n", - "I_0 = IdentityOperator(V0_h.vector_space)\n", + "I_0 = IdentityOperator(V0_h.coeff_space)\n", "P_Gamma = I_0 - P_0" ] }, diff --git a/exercises/poisson_2d_I_naturalBC.ipynb b/exercises/poisson_2d_I_naturalBC.ipynb index da5fabe..58e6866 100644 --- a/exercises/poisson_2d_I_naturalBC.ipynb +++ b/exercises/poisson_2d_I_naturalBC.ipynb @@ -1,12 +1,5 @@ { "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": {}, diff --git a/exercises/time_harmonic_maxwell_2d_I.ipynb b/exercises/time_harmonic_maxwell_2d_I.ipynb index 5d2c5a6..98be3e1 100644 --- a/exercises/time_harmonic_maxwell_2d_I.ipynb +++ b/exercises/time_harmonic_maxwell_2d_I.ipynb @@ -145,9 +145,9 @@ "from utils import HcurlBoundaryProjector2D\n", "from psydac.linalg.basic import IdentityOperator\n", "\n", - "P_Hcurl = HcurlBoundaryProjector2D(V, V_h.vector_space)\n", + "P_Hcurl = HcurlBoundaryProjector2D(V, V_h.coeff_space)\n", "\n", - "I1 = IdentityOperator(V_h.vector_space)\n", + "I1 = IdentityOperator(V_h.coeff_space)\n", "P_Gamma = I1 - P_Hcurl" ] }, diff --git a/exercises/time_harmonic_maxwell_2d_II.ipynb b/exercises/time_harmonic_maxwell_2d_II.ipynb index 2a84277..abf33ef 100644 --- a/exercises/time_harmonic_maxwell_2d_II.ipynb +++ b/exercises/time_harmonic_maxwell_2d_II.ipynb @@ -156,9 +156,9 @@ "from utils import HcurlBoundaryProjector2D\n", "from psydac.linalg.basic import IdentityOperator\n", "\n", - "P_Hcurl = HcurlBoundaryProjector2D(V1, V1_h.vector_space)\n", + "P_Hcurl = HcurlBoundaryProjector2D(V1, V1_h.coeff_space)\n", "\n", - "I1 = IdentityOperator(V1_h.vector_space)\n", + "I1 = IdentityOperator(V1_h.coeff_space)\n", "P_Gamma = I1 - P_Hcurl" ] }, From 15c5a39bed4a46e3cd44919f878b1224ddc4c8d5 Mon Sep 17 00:00:00 2001 From: jowezarek Date: Fri, 4 Apr 2025 10:45:57 +0200 Subject: [PATCH 16/16] minor tidy up --- exercises/electro_static_2d.ipynb | 64 ++++++++++++++++--------------- 1 file changed, 34 insertions(+), 30 deletions(-) diff --git a/exercises/electro_static_2d.ipynb b/exercises/electro_static_2d.ipynb index 399edce..0a7b1d5 100644 --- a/exercises/electro_static_2d.ipynb +++ b/exercises/electro_static_2d.ipynb @@ -143,21 +143,24 @@ "#P_H1_Gamma = I0 - P_H1\n", "#P_Hcurl_Gamma = I1 - P_Hcurl\n", "\n", - "# Modified Operators for Projection Method (using P_H1, P_Hcurl, P_H1_Gamma, P_Hcurl_Gamma)\n", - "# <-\n", - "# ...\n", - "# <-\n", - "\n", - "# Inverse M0 mass matrix\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", "\n", - "# System Matrix A\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\n", + "# rhs vector f and f_bc\n", "#rho_coeffs = P0(rho_lambdified).coeffs\n", - "#f = # <-" + "#f = # <-\n", + "#f_bc = P_Hcurl @ f" ] }, { @@ -171,13 +174,15 @@ "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_inv = inverse(A, 'pcg', pc=M1_0_inv, tol=tol, maxiter=maxiter)\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_inv @ f\n", + "#E_h = A_bc_inv @ f_bc\n", "t1 = time.time()" ] }, @@ -226,7 +231,7 @@ "# 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_inv.get_info() )\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 ) )" @@ -256,26 +261,25 @@ "#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", - "if False:\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", - " )" + "#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 + "nbformat_minor": 4 }