-
Notifications
You must be signed in to change notification settings - Fork 3
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
1 parent
1adfe25
commit a992d1a
Showing
14 changed files
with
1,133 additions
and
38 deletions.
There are no files selected for viewing
Binary file added
BIN
+5.9 KB
_images/04da5cd95c68db0373a0680a8ff38a87d14cc0b24722dda60ce7c4a5dc37b760.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -1,2 +1,10 @@ | ||
| (examples)= | ||
| ## Examples | ||
| # Examples | ||
|
|
||
| ```{contents} | ||
| :local: | ||
| ``` | ||
|
|
||
|
|
||
| ```{tableofcontents} | ||
| ``` |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,310 @@ | ||
| { | ||
| "cells": [ | ||
| { | ||
| "cell_type": "markdown", | ||
| "id": "98d572eb", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "# A plane wave on a square\n", | ||
| "\n", | ||
| "We solve the two-dimensional wave equation to find $H: [0,T]\\to H^1(\\Omega)$ and the vector field $E: [0,T]\\to H(\\mathrm{div})$ are\n", | ||
| "\n", | ||
| "$$\n", | ||
| "\\begin{aligned}\n", | ||
| "\\partial_t E &= -\\nabla H,&t\\in(0,T),x\\in\\Omega\\\\\n", | ||
| "\\partial_t H&= -\\mathrm{div} E + f,&t\\in(0,T),x\\in\\Omega\\\\\n", | ||
| "H(0,x) &= \\exp(-400(y-1/2)^2),&t\\in(0,T),x\\in\\Omega\\\\\n", | ||
| "E(0,x) &= 0,&x\\in\\Omega\\\\\n", | ||
| "H(t,x) &= 0,&t\\in(0,T),x\\in\\partial\\Omega\n", | ||
| "\\end{aligned}\n", | ||
| "$$\n", | ||
| "for $\\Omega=(0,1)^2$\n" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": 1, | ||
| "id": "ef5e1523", | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "from ngsolve import *\n", | ||
| "import dualcellspaces as dcs\n", | ||
| "from time import time\n", | ||
| "from ngsolve.webgui import Draw" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "id": "9f4aa8f6", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "After the necessary imports we define some parameters and the mesh" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": 12, | ||
| "id": "4b18a2aa", | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "maxh = 0.1\n", | ||
| "tend = 3\n", | ||
| "order = 2\n", | ||
| "\n", | ||
| "H0 = CF(exp(-20**2*((y-1/2)**2)))\n", | ||
| "E0 = CF((0,0))\n", | ||
| "\n", | ||
| "mesh = Mesh(unit_square.GenerateMesh(maxh=maxh))" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "id": "d11b5b72", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "We define the spaces" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": 13, | ||
| "id": "d1052349", | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "fesH = dcs.H1DualCells(mesh, order=order)\n", | ||
| "fesE = dcs.HDivPrimalCells(mesh, order=order)" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "id": "5e6daede", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "To define the (Mass) bilinear forms we need special integration rules:" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": 14, | ||
| "id": "448f75fa", | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "E, dE = fesE.TnT()\n", | ||
| "H, dH = fesH.TnT()\n", | ||
| "\n", | ||
| "dxH = dx(intrules=fesH.GetIntegrationRules())\n", | ||
| "dSw = dx(element_boundary=True,intrules=dcs.GetIntegrationRules(2*order+6))\n", | ||
| "dxw = dx(intrules=dcs.GetIntegrationRules(2*order+6))\n", | ||
| "\n", | ||
| "\n", | ||
| "massE = fesE.Mass(Id(2))\n", | ||
| "massH = fesH.Mass(1)\n", | ||
| "massinvE = massE.Inverse()\n", | ||
| "massinvH = massH.Inverse()\n", | ||
| "\n", | ||
| "normal = specialcf.normal(2)\n", | ||
| "\n", | ||
| "Grad = BilinearForm(-H*div(dE)*dxw+H*dE*normal*dSw, geom_free=True).Assemble().mat\n", | ||
| "\n", | ||
| "lffH = LinearForm(dH*H0*dxH).Assemble()" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "id": "513db898", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "The maximal admissible time step may be estimated using a simple power iteration" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "id": "d921ecd3", | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "def estimate_tau(mat, maxsteps = 1000, tol = 1e-4): \n", | ||
| " vec = mat.CreateColVector()\n", | ||
| " vec.SetRandom()\n", | ||
| " tmp = vec.CreateVector()\n", | ||
| " lam = 0\n", | ||
| " for i in range(maxsteps):\n", | ||
| " #print(i,end='\\r')\n", | ||
| " tmp.data = mat * vec\n", | ||
| " \n", | ||
| " lamnew = InnerProduct(tmp,vec)\n", | ||
| " tau = 2/sqrt(lamnew)\n", | ||
| " #res=(lamnew*vec-tmp).Norm()\n", | ||
| " tmp *= 1/tmp.Norm()\n", | ||
| " #print(lamnew)\n", | ||
| " diff = (tmp-vec).Norm()\n", | ||
| " if diff<tol: return tau\n", | ||
| " vec.data = tmp\n", | ||
| " lam = lamnew\n", | ||
| " print(\"did not converge, last diff = \",diff)\n", | ||
| " return tau" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": 5, | ||
| "id": "b519eeab", | ||
| "metadata": {}, | ||
| "outputs": [ | ||
| { | ||
| "name": "stdout", | ||
| "output_type": "stream", | ||
| "text": [ | ||
| "estimated timestep tau: 4.678488e-03\n" | ||
| ] | ||
| } | ||
| ], | ||
| "source": [ | ||
| "tau = estimate_tau([email protected]@massinvE@Grad)\n", | ||
| "\n", | ||
| "print(\"estimated timestep tau: {:e}\".format(tau))\n", | ||
| "tau*=0.9" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "id": "1c211d43", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "It remains to set the initial conditions..." | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "id": "5b30c8aa", | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "t = 0.\n", | ||
| "i = 0\n", | ||
| "drawevery = 10\n", | ||
| "\n", | ||
| "gfE = GridFunction(fesE)\n", | ||
| "gfH = GridFunction(fesH)\n", | ||
| "\n", | ||
| "gfH_history = GridFunction(fesH,multidim=0)\n", | ||
| "\n", | ||
| "gfH.vec.data = massinvH*lffH.vec\n", | ||
| "gfE.vec.data[:] = 0.\n", | ||
| "\n", | ||
| "#scene = Draw(gfH,mesh,intpoints=dcs.GetWebGuiPoints(2),order=2,autoscale=False,min=0,max=1)\n", | ||
| "\n", | ||
| "gfE.vec.data = tau/2*massinvE@Grad*gfH.vec" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "id": "adad4709", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "... and start the time loop" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": 7, | ||
| "id": "06f43b80", | ||
| "metadata": {}, | ||
| "outputs": [ | ||
| { | ||
| "name": "stdout", | ||
| "output_type": "stream", | ||
| "text": [ | ||
| " time = 2.989553877139938\t step = 710\t energy = 0.06208402112680203\t current dofs/s = 7.691518e+0666\n", | ||
| " timesteps: 712.4808876291984\t dofs: 14452\t dofs per second: 7.972132e+06\n" | ||
| ] | ||
| } | ||
| ], | ||
| "source": [ | ||
| "\n", | ||
| "now = time()\n", | ||
| "nowstart = now\n", | ||
| "\n", | ||
| "times = []\n", | ||
| "energies = []\n", | ||
| "tmpH = gfH.vec.CreateVector()\n", | ||
| "tmpE = gfE.vec.CreateVector()\n", | ||
| "with TaskManager():\n", | ||
| " while t<tend:\n", | ||
| " if i%drawevery == 0:\n", | ||
| " gfH_history.AddMultiDimComponent(gfH.vec)\n", | ||
| " #scene.Redraw()\n", | ||
| " times.append(t)\n", | ||
| " tmpH.data = massH * gfH.vec\n", | ||
| " tmpE.data = massE * gfE.vec\n", | ||
| " energies.append(InnerProduct(gfE.vec,tmpE)+InnerProduct(gfH.vec,tmpH))\n", | ||
| " #print(\"\\r time = {}\\t step = {}\\t energy = {}\\t current dofs/s = {:e}\".format(t,i,energies[-1],(fesE.ndof+fesH.ndof)*drawevery/(time()-now)),end=\"\")\n", | ||
| " now = time()\n", | ||
| " i=i+1\n", | ||
| " t+=tau\n", | ||
| " gfH.vec.data += -tau*[email protected]*gfE.vec\n", | ||
| " gfE.vec.data += tau*massinvE@Grad*gfH.vec\n", | ||
| "\n", | ||
| "comptime = time()-nowstart\n", | ||
| "print(\"\\n timesteps: {}\\t dofs: {}\\t dofs per second: {:e}\".format(tend/tau, (fesE.ndof+fesH.ndof),(fesE.ndof+fesH.ndof)*tend/tau/comptime))" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "id": "e85c89eb", | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "scene = Draw(gfH_history,mesh,intpoints=dcs.GetWebGuiPoints(2),order=2,autoscale=False,min=0,max=1,animate=True)" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "id": "1dbae73d", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "We observe preservation of a modified (discrete) energy:" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "id": "bb7a7edd", | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "import matplotlib.pyplot as pl\n", | ||
| "pl.plot(times,energies)\n", | ||
| "pl.ylim((0,0.1))" | ||
| ] | ||
| } | ||
| ], | ||
| "metadata": { | ||
| "kernelspec": { | ||
| "display_name": "Python 3 (ipykernel)", | ||
| "language": "python", | ||
| "name": "python3" | ||
| }, | ||
| "language_info": { | ||
| "codemirror_mode": { | ||
| "name": "ipython", | ||
| "version": 3 | ||
| }, | ||
| "file_extension": ".py", | ||
| "mimetype": "text/x-python", | ||
| "name": "python", | ||
| "nbconvert_exporter": "python", | ||
| "pygments_lexer": "ipython3", | ||
| "version": "3.10.12" | ||
| } | ||
| }, | ||
| "nbformat": 4, | ||
| "nbformat_minor": 5 | ||
| } |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Oops, something went wrong.