diff --git a/_toc.yml b/_toc.yml index b98b174..0804ebd 100644 --- a/_toc.yml +++ b/_toc.yml @@ -33,6 +33,7 @@ parts: - file: chapter1/space - file: chapter1/rules - file: chapter1/poisson + - file: chapter1/boundary-conditions - caption: Linear Problems chapters: - file: chapter2/poisson diff --git a/chapter1/boundary-conditions.md b/chapter1/boundary-conditions.md new file mode 100644 index 0000000..6764b3d --- /dev/null +++ b/chapter1/boundary-conditions.md @@ -0,0 +1,64 @@ +# Boundary Conditions +*Author: Ahmed Ratnani* + +SymPDE & Psydac allows you to use both strong and weak boundary conditions. +We start first by explaining how to identify a boundary in **NCube** domains, such as **Line**, **Square** and a **Cube**. + +## Identifying a boundary in **NCube** domains + +Since we're dealing with **NCube** domains, the best way to identify a boundary is to use the couple **(axis, extremity)**. The axis is usually defined as perpendicular to the boundary, while the **extremity** gives the lower or upper bound of the interval associated to the **axis**. + +Let's see it through the following examples; + +For a **Line**, the following figure gives the value to access a boundary. + +![png](images/boundary-conditions/line-boundaries.png) + +For a **Square**, we have, + +![png](images/boundary-conditions/square-boundaries.png) + +For a **Cube**, it is similar and would not be helpful to visualize it here. + +## Essential boundary conditions + +Essential boundary conditions are usually treated in two ways, either in the strong or weak form. The latter can be achieved through the use of Nitsche's method. You can check the associated examples in the sequel. + +For the essential boundary conditions, usually, one needs one of the following expressions; + +$$ +\begin{align} +u &= 0, \quad \text{(homogeneous case)}, \\ +u &= f, \quad \text{(inhomogeneous case)}, \\ +\mathbf{v} &= 0, \quad \text{(homogeneous case)}, \\ +\mathbf{v} &= \mathbf{g}, \quad \text{(homogeneous case)}, +\end{align} +$$ + +where $u$ and $f$ are scalar functions while $\mathbf{v}$ and $\mathbf{g}$ denote vector functions. + + +Let `Gamma` denotes the boundary $\Gamma$. +The normal derivative is an available operator in SymPDE, it is provided as `Dn` operator. Alternatively, you can also define it manually, + +```python + nn = NormalVector('nn') + dn = lambda a: dot(grad(a), nn) +``` + + +we have the following equivalences, + +| Mathematical Expression | SymPDE Expression | Example | +| --------------------------------------------------------- | ------------------------------ | ---------------------------------- | +| $u = 0$ on $\Gamma$ | `EssentialBC(u, 0, Gamma)` | | +| $u = f$ on $\Gamma$ | `EssentialBC(u, f, Gamma)` | | +| $\partial_n u = 0$ on $\Gamma$ | `EssentialBC(dn(u), 0, Gamma)` | | +| $\partial_n u = f$ on $\Gamma$ | `EssentialBC(dn(u), f, Gamma)` | | +| $\mathbf{v} = 0$ on $\Gamma$ | `EssentialBC(v, 0, Gamma)` | | +| $\mathbf{v} = \mathbf{g}$ on $\Gamma$ | `EssentialBC(v, g, Gamma)` | | +| $\mathbf{v} \cdot \mathbf{n} = 0$ on $\Gamma$ | ? | | +| $\mathbf{v} \times \mathbf{n} = 0$ on $\Gamma$ | ? | | +| $\mathbf{v} \cdot \mathbf{n} = f,\mathbf{g}$ on $\Gamma$ | ? | | +| $\mathbf{v} \times \mathbf{n} = f,\mathbf{g}$ on $\Gamma$ | ? | | + diff --git a/chapter1/images/boundary-conditions/line-boundaries.png b/chapter1/images/boundary-conditions/line-boundaries.png new file mode 100644 index 0000000..b5e38ba Binary files /dev/null and b/chapter1/images/boundary-conditions/line-boundaries.png differ diff --git a/chapter1/images/boundary-conditions/square-boundaries.png b/chapter1/images/boundary-conditions/square-boundaries.png new file mode 100644 index 0000000..818c117 Binary files /dev/null and b/chapter1/images/boundary-conditions/square-boundaries.png differ diff --git a/chapter2/poisson-nitsche.ipynb b/chapter2/poisson-nitsche.ipynb index 052e688..68d5626 100644 --- a/chapter2/poisson-nitsche.ipynb +++ b/chapter2/poisson-nitsche.ipynb @@ -18,8 +18,8 @@ "\n", "$$\n", "\\begin{align}\n", - " - \\nabla^2 u = f \\quad \\text{in~$\\Omega$}, \\quad \\quad \n", - " u = g \\quad \\text{on~$\\partial \\Omega$}.\n", + " - \\nabla^2 u = f \\quad \\text{in $\\Omega$}, \\quad \\quad \n", + " u = g \\quad \\text{on $\\partial \\Omega$}.\n", "\\end{align}\n", "$$" ] @@ -54,12 +54,12 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 18, "id": "ecd1330a", "metadata": {}, "outputs": [], "source": [ - "from sympde.expr import BilinearForm, LinearForm, integral\n", + "from sympde.expr import BilinearForm, LinearForm, integral, Norm\n", "from sympde.expr import find, EssentialBC\n", "from sympde.topology import ScalarFunctionSpace, Square, element_of\n", "from sympde.calculus import grad, dot, Dn\n", @@ -80,16 +80,19 @@ "u,v = [element_of(V, name=i) for i in ['u', 'v']]\n", "\n", "# bilinear form\n", - "a = BilinearForm((u,v), integral(domain, dot(grad(v), grad(u)) + kappa*u*v - u*Dn(v) - v*Dn(u)))\n", + "a = BilinearForm((u,v), integral(domain, dot(grad(v), grad(u))))\n", + "an = BilinearForm((u,v), integral(domain.boundary, kappa*u*v - u*Dn(v) - v*Dn(u)))\n", "\n", "# linear form\n", - "g = sin(pi*x)*sin(pi*y)\n", - "f = 2*pi**2*sin(pi*x)*sin(pi*y)\n", + "ue = sin(pi*x)*sin(pi*y)\n", + "g = ue\n", + "f = 2*pi**2*sin(pi*x)*sin(pi*y)\n", "\n", - "l = LinearForm(v, integral(domain, f*v + kappa*g*v - g*Dn(v)))\n", + "l = LinearForm(v, integral(domain, f*v))\n", + "ln = LinearForm(v, integral(domain.boundary, kappa*g*v - g*Dn(v)))\n", "\n", "# Variational problem\n", - "equation = find(u, forall=v, lhs=a(u, v), rhs=l(v))" + "equation = find(u, forall=v, lhs=a(u,v) + an(u,v), rhs=l(v) + ln(v))" ] }, { @@ -102,7 +105,7 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 19, "id": "3d86674b", "metadata": {}, "outputs": [], @@ -113,7 +116,7 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 20, "id": "94c83267", "metadata": {}, "outputs": [], @@ -125,7 +128,7 @@ "Vh = discretize(V, domain_h, degree=degree)\n", "\n", "# Discretize equation\n", - "#equation_h = discretize(equation, domain_h, [Vh, Vh])" + "equation_h = discretize(equation, domain_h, [Vh, Vh])" ] }, { @@ -138,53 +141,69 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 21, "id": "8f1b2d68", "metadata": {}, "outputs": [], "source": [ - "#uh = equation_h.solve()" + "equation_h.set_solver('gmres', info=False, tol=1e-8)" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 22, "id": "e96f321c", "metadata": {}, "outputs": [], - "source": [] + "source": [ + "uh = equation_h.solve(kappa=1e3)" + ] }, { - "cell_type": "code", - "execution_count": null, - "id": "c7ae0b9e", + "cell_type": "markdown", + "id": "d87fd79f", "metadata": {}, - "outputs": [], - "source": [] + "source": [ + "## Computing the error norm" + ] }, { - "cell_type": "code", - "execution_count": null, - "id": "9ad386a6", + "cell_type": "markdown", + "id": "08ee5cac", "metadata": {}, - "outputs": [], - "source": [] + "source": [ + "### Computing the $L^2$ norm" + ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 23, "id": "068bdb95", "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "aa88a918", - "metadata": {}, - "outputs": [], - "source": [] + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0.00021402394834116976\n" + ] + } + ], + "source": [ + "u = element_of(V, name='u')\n", + "\n", + "# create the formal Norm object\n", + "l2norm = Norm(u - ue, 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)" + ] } ], "metadata": { diff --git a/chapter2/poisson.ipynb b/chapter2/poisson.ipynb index de5c678..0e26eb1 100644 --- a/chapter2/poisson.ipynb +++ b/chapter2/poisson.ipynb @@ -12,9 +12,10 @@ "\n", "$$\n", "\\begin{align}\n", - " - \\nabla^2 u = f \\quad \\text{in $\\Omega$}, \\quad \\quad \n", - " u = 0 \\quad \\text{on $\\Gamma_D$}, \\quad \\quad \n", - " \\partial_n u = g \\quad \\text{on $\\Gamma_N$}.\n", + " - \\nabla^2 u = f \\quad &\\text{in $\\Omega$}, \\\\ \n", + " u = 0 \\quad &\\text{on $\\Gamma_0$}, \\\\\n", + " u = g_i \\quad &\\text{on $\\Gamma_I$}, \\\\\n", + " \\partial_n u = g_n \\quad &\\text{on $\\Gamma_N := \\partial \\Omega \\setminus \\left( \\Gamma_0 \\cup \\Gamma_I \\right)$}.\n", "\\end{align}\n", "$$\n", "\n", @@ -32,13 +33,7 @@ "\n", "- $V \\subset H^1(\\Omega)$, \n", "- $a(u,v) := \\int_{\\Omega} \\nabla u \\cdot \\nabla v ~ d\\Omega$,\n", - "- $l(v) := \\int_{\\Omega} f v ~ d\\Omega + \\int_{\\Gamma_N} g v ~ d\\Gamma$.\n", - "\n", - "## Examples\n", - "\n", - "- [Poisson on a square domain with Homogeneous Boundary Conditions](https://github.com/pyccel/sympde/blob/devel-documentation/docs/examples/2d_poisson_dir0.py)\n", - "- [Poisson on a square domain with Homogeneous and Neumann Boundary Conditions](https://github.com/pyccel/sympde/blob/devel-documentation/docs/examples/2d_poisson_dir0_neu.py)\n", - "- [Poisson on a square domain with Homogeneous, Inhomogeneous and Neumann Boundary Conditions](https://github.com/pyccel/sympde/blob/devel-documentation/docs/examples/2d_poisson_dir_neu.py)" + "- $l(v) := \\int_{\\Omega} f v ~ d\\Omega + \\int_{\\Gamma_N} g_n v ~ d\\Gamma$." ] }, { @@ -51,19 +46,26 @@ }, { "cell_type": "code", - "execution_count": 27, + "execution_count": 8, "id": "d742586c", "metadata": {}, "outputs": [], "source": [ - "from sympde.expr import BilinearForm, LinearForm, integral\n", + "from sympde.expr import BilinearForm, LinearForm, integral\n", "from sympde.expr import find, EssentialBC, Norm, SemiNorm\n", "from sympde.topology import ScalarFunctionSpace, Square, element_of\n", - "from sympde.calculus import grad, dot\n", + "from sympde.calculus import grad, dot, laplace\n", + "from sympde.topology import NormalVector, Union\n", "\n", - "from sympy import pi, sin\n", + "from sympy import pi, sin\n", + "\n", + "from psydac.api.discretization import discretize\n", "\n", "domain = Square()\n", + "Gamma_0 = domain.get_boundary(axis=0, ext=-1)\n", + "Gamma_i = domain.get_boundary(axis=0, ext=1)\n", + "Gamma_n = domain.boundary.complement(Union(Gamma_0, Gamma_i))\n", + "nn = NormalVector('nn')\n", "\n", "V = ScalarFunctionSpace('V', domain)\n", "\n", @@ -74,34 +76,25 @@ "# bilinear form\n", "a = BilinearForm((u,v), integral(domain , dot(grad(v), grad(u))))\n", "\n", + "# exact solution\n", + "ue = sin(pi*x) * (1+y*sin(pi*y/3))**2\n", + "L = lambda w: - laplace(w)\n", + "f = L(ue)\n", + "gi = ue\n", + "gn = ue\n", + "\n", "# linear form\n", - "f = 2*pi**2*sin(pi*x)*sin(pi*y)\n", "l = LinearForm(v, integral(domain, f*v))\n", "\n", + "# Boundary term for the Neumann BC\n", + "ln = LinearForm(v, integral(Gamma_n, v * dot(grad(gn), nn)))\n", + "\n", "# Dirichlet boundary conditions\n", - "bc = [EssentialBC(u, 0, domain.boundary)]\n", + "bc = [EssentialBC(u, 0, Gamma_0)]\n", + "bc += [EssentialBC(u, gi, Gamma_i)]\n", "\n", "# Variational problem\n", - "equation = find(u, forall=v, lhs=a(u, v), rhs=l(v), bc=bc)" - ] - }, - { - "cell_type": "markdown", - "id": "62ac1fd4", - "metadata": {}, - "source": [ - "\n", - "This very simple Python code reflects well the abstract mathematical framework needed for variational formulations.\n", - "The structure of the code is as follows,\n", - "\n", - "1. Create a domain.\n", - "2. Create a space of *scalar* functions over the domain.\n", - "3. Create elements from this function space. These elements will denote the test and trial functions.\n", - "4. Create the Bilinear and Linear forms, $a$ and $l$ respectively.\n", - "5. Create Essential Boundary Conditions.\n", - "6. Create the variational problem.\n", - "\n", - "Most of the time, you will need to follow the same steps, with some minor variants depending on the problem you're considering." + "equation = find(u, forall=v, lhs=a(u, v), rhs=l(v)+ln(v), bc=bc)" ] }, { @@ -112,27 +105,9 @@ "## Discretization" ] }, - { - "cell_type": "markdown", - "id": "51095918", - "metadata": {}, - "source": [ - "We shall need the **discretize** function from **PsyDAC**." - ] - }, - { - "cell_type": "code", - "execution_count": 21, - "id": "a2a0a2a1", - "metadata": {}, - "outputs": [], - "source": [ - "from psydac.api.discretization import discretize" - ] - }, { "cell_type": "code", - "execution_count": 22, + "execution_count": 23, "id": "00e54163", "metadata": {}, "outputs": [], @@ -143,7 +118,7 @@ }, { "cell_type": "code", - "execution_count": 23, + "execution_count": 24, "id": "5999c62b", "metadata": {}, "outputs": [], @@ -168,7 +143,17 @@ }, { "cell_type": "code", - "execution_count": 24, + "execution_count": 25, + "id": "004dfdb3", + "metadata": {}, + "outputs": [], + "source": [ + "equation_h.set_solver('gmres', info=False, tol=1e-8)" + ] + }, + { + "cell_type": "code", + "execution_count": 26, "id": "541192ee", "metadata": {}, "outputs": [], @@ -202,7 +187,7 @@ }, { "cell_type": "code", - "execution_count": 26, + "execution_count": 27, "id": "5925c6cd", "metadata": {}, "outputs": [ @@ -210,13 +195,11 @@ "name": "stdout", "output_type": "stream", "text": [ - "9.49756716724664e-07\n" + "0.00042499840633644024\n" ] } ], "source": [ - "ue = sin(pi*x)*sin(pi*y)\n", - "\n", "u = element_of(V, name='u')\n", "\n", "# create the formal Norm object\n", @@ -242,7 +225,7 @@ }, { "cell_type": "code", - "execution_count": 30, + "execution_count": 28, "id": "e5c1a8b8", "metadata": {}, "outputs": [ @@ -250,7 +233,7 @@ "name": "stdout", "output_type": "stream", "text": [ - "9.768702410721585e-05\n" + "0.025283770887649267\n" ] } ], @@ -286,7 +269,7 @@ "name": "stdout", "output_type": "stream", "text": [ - "9.769164089397408e-05\n" + "0.025287342563909892\n" ] } ], @@ -303,14 +286,6 @@ "# print the result\n", "print(h1_error)" ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "9250897b", - "metadata": {}, - "outputs": [], - "source": [] } ], "metadata": { diff --git a/chapter4/subdomains.md b/chapter4/subdomains.md index ef83d5e..51cb34d 100644 --- a/chapter4/subdomains.md +++ b/chapter4/subdomains.md @@ -1,4 +1,5 @@ # Subdomains +*Author: Ahmed Ratnani* In this section, we consider a domain $\Omega$ which a union of multiple subdomain, *i.e.*