Skip to content

Commit

Permalink
update notebooks
Browse files Browse the repository at this point in the history
  • Loading branch information
@PaulSt
PaulSt committed Jun 11, 2024
1 parent 8ffda7f commit c0ab962
Show file tree
Hide file tree
Showing 5 changed files with 98 additions and 99 deletions.
138 changes: 82 additions & 56 deletions docs/notebooks/embTrefftz-helm.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,8 +55,8 @@
"\\newcommand{\\inner}[1]{\\langle #1 \\rangle}\n",
"\\begin{align*}\n",
" \\begin{cases}\n",
" -\\Delta u - \\omega^2 u= 0 &\\text{ in } \\dom, \\\\\n",
" \\frac{\\partial u}{\\partial \\nx} + i u = g &\\text{ on } \\partial \\dom.\n",
" -\\Delta u - \\omega^2 u= f &\\text{ in } \\dom, \\\\\n",
" \\frac{\\partial u}{\\partial \\nx} - i\\omega u = g &\\text{ on } \\partial \\dom.\n",
" \\end{cases}\n",
"\\end{align*}\n",
"$$\n",
Expand Down Expand Up @@ -90,15 +90,15 @@
"metadata": {},
"outputs": [],
"source": [
"def TrefftzHelmholtzEmb(fes):\n",
"def TrefftzHelmholtzEmb(fes,omega=1):\n",
" mesh = fes.mesh\n",
" k = fes.globalorder\n",
" u = fes.TrialFunction()\n",
" \n",
" Lap = lambda u : sum(Trace(u.Operator('hesse')))\n",
" fes2 = L2(mesh, order=order-2, dgjumps=True,complex=True)\n",
" v = fes2.TestFunction()\n",
" op = -u*v*dx - Lap(u)*v*dx\n",
" op = -omega**2*u*v*dx - Lap(u)*v*dx\n",
" PP = TrefftzEmbedding(op,fes,test_fes=fes2)\n",
" return PP"
]
Expand Down Expand Up @@ -163,19 +163,6 @@
"Draw (gfshow, mesh, interpolate_multidim=False, animate=False, autoscale=False, min=0.8,max=1, settings={\"subdivision\":20})"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ad81be3d",
"metadata": {},
"outputs": [],
"source": [
"mesh = Mesh(unit_square.GenerateMesh(maxh=.3))\n",
"order = 4\n",
"fes = L2(mesh, order=order, complex=True, dgjumps=True)\n",
"u,v = fes.TnT()"
]
},
{
"cell_type": "code",
"execution_count": null,
Expand All @@ -187,7 +174,7 @@
"n = specialcf.normal(2)\n",
"exact = exp(1j*sqrt(0.5)*(x+y))\n",
"gradexact = CoefficientFunction((sqrt(0.5)*1j*exact, sqrt(0.5)*1j*exact))\n",
"bndc = gradexact*n + 1j*omega*exact\n",
"bndc = gradexact*n - 1j*omega*exact\n",
"eps = 10**-7"
]
},
Expand All @@ -202,9 +189,9 @@
"\\begin{align}\n",
" a_h(u,v) &= \\sum_{K\\in\\Th}\\int_K \\nabla u\\nabla v-\\omega^2 uv\\ dV\n",
" -\\int_{\\Fh^\\text{int}}\\left(\\avg{\\nabla u}\\jump{v}+\\jump{u} \\avg{\\overline{\\nabla v}} \\right) dS \\nonumber \\\\ \n",
" &\\qquad+\\int_{\\Fh^\\text{int}} \\left( i\\alpha \\omega\\jump{u}\\jump{\\overline{v}} - \\frac{\\beta}{i\\omega}\\jump{\\nabla u}\\jump{\\overline{\\nabla v}} \\right) dS -\\int_{\\Fh^\\text{bnd}}\\delta\\left(\\nx\\cdot\\nabla u \\overline{v}+u \\overline{\\nx\\cdot\\nabla v}\\right) dS\\\\ \\nonumber\n",
" &\\qquad+\\int_{\\Fh^\\text{bnd}} \\left( i(1-\\delta)\\omega{u}{\\overline{v}} - \\frac{\\delta}{i\\omega}{\\nabla u}{\\overline{\\nabla v}} \\right) dS \\\\ \n",
" \\ell(v) &= \\int_{\\Fh^\\text{bnd}}\\left( (1-\\delta)g\\overline{v} - \\frac{\\delta}{i\\omega}g\\overline{\\nx\\cdot\\nabla v}\\right) dS\n",
" &\\qquad+\\int_{\\Fh^\\text{int}} \\left( i\\alpha \\omega\\jump{u}\\jump{\\overline{v}} - i\\frac{\\beta}{\\omega}\\jump{\\nabla u}\\jump{\\overline{\\nabla v}} \\right) dS -\\int_{\\Fh^\\text{bnd}}\\delta\\left(\\nx\\cdot\\nabla u \\overline{v}+u \\overline{\\nx\\cdot\\nabla v}\\right) dS\\\\ \\nonumber\n",
" &\\qquad-\\int_{\\Fh^\\text{bnd}} \\left( i(1-\\delta)\\omega{u}{\\overline{v}} + i\\frac{\\delta}{\\omega}{\\nabla u\\cdot\\nx}{\\overline{\\nabla v}\\cdot\\nx} \\right) dS \\\\ \n",
" \\ell(v) &= \\int_\\Omega fv\\ dV + \\int_{\\Fh^\\text{bnd}}\\left( (1-\\delta)g\\overline{v} - i\\frac{\\delta}{\\omega}g\\overline{\\nx\\cdot\\nabla v}\\right) dS\n",
"\\end{align}\n",
"$$"
]
Expand All @@ -216,37 +203,43 @@
"metadata": {},
"outputs": [],
"source": [
"h = specialcf.mesh_size\n",
"alpha = 1/(omega*h)\n",
"beta = omega*h\n",
"delta = omega*h\n",
"\n",
"jump_u = (u-u.Other())*n\n",
"jump_v = (v-v.Other())*n\n",
"jump_du = (grad(u)-grad(u.Other()))*n\n",
"jump_dv = (grad(v)-grad(v.Other()))*n\n",
"mean_u = 0.5 * ((u)+(u.Other()))\n",
"mean_du = 0.5 * (grad(u)+grad(u.Other()))\n",
"mean_dv = 0.5 * (grad(v)+grad(v.Other()))\n",
"\n",
"a = BilinearForm(fes)\n",
"a += grad(u)*(grad(v))*dx - omega**2*u*(v)*dx\n",
"\n",
"a += -(jump_u*(mean_dv)+mean_du*(jump_v)) * dx(skeleton=True)\n",
"a += -1/(omega*1j)*beta*(jump_du*(jump_dv)) * dx(skeleton=True)\n",
"a += omega*1j*alpha*jump_u*(jump_v) * dx(skeleton=True)\n",
"\n",
"a += -delta*(u*(grad(v))*n+grad(u)*n*(v)) * ds(skeleton=True)\n",
"a += -1/(omega*1j)*delta*(grad(u)*n)*((grad(v))*n) * ds(skeleton=True)\n",
"a += omega*1j*(1-delta)*u*(v) * ds(skeleton=True)\n",
"\n",
"f = LinearForm(fes)\n",
"f += -1/(omega*1j)*delta*bndc*(grad(v))*n*ds(skeleton=True)\n",
"f += (1-delta)*bndc*(v)*ds(skeleton=True)\n",
"\n",
"with TaskManager():\n",
" a.Assemble()\n",
" f.Assemble()"
"def dghelmholtz(fes,omega,rhs=0,bndc=0,bnd=\".*\"):\n",
" u,v = fes.TnT()\n",
" h = specialcf.mesh_size\n",
" order = fes.globalorder\n",
" alpha = order**2/h\n",
" beta = h/order\n",
" delta = 0.1*omega*h/order\n",
" IP = lambda u,v: InnerProduct(u,v)\n",
" \n",
" jump_u = (u-u.Other())*n\n",
" jump_v = (v-v.Other())*n\n",
" jump_du = (grad(u)-grad(u.Other()))*n\n",
" jump_dv = (grad(v)-grad(v.Other()))*n\n",
" mean_u = 0.5 * (u+u.Other())\n",
" mean_du = 0.5 * (grad(u)+grad(u.Other()))\n",
" mean_dv = 0.5 * (grad(v)+grad(v.Other()))\n",
" \n",
" a = BilinearForm(fes)\n",
" a += grad(u)*grad(v)*dx - omega**2*u*v*dx\n",
" \n",
" a += -(jump_u*mean_dv+mean_du*jump_v) * dx(skeleton=True)\n",
" a += -1j*omega*alpha*jump_u*jump_v * dx(skeleton=True)\n",
" a += -1j*beta/(omega)*(jump_du*jump_dv) * dx(skeleton=True)\n",
" \n",
" a += -delta*(u*grad(v)*n+grad(u)*n*v) * ds(bnd,skeleton=True)\n",
" a += -1j*(1-delta)*omega*u*v * ds(bnd,skeleton=True)\n",
" a += -1j*delta/(omega)*(grad(u)*n)*(grad(v)*n) * ds(bnd,skeleton=True)\n",
" \n",
" f = LinearForm(fes)\n",
" f += rhs * v * dx\n",
" f += (1-delta)*bndc*v * ds(bnd,skeleton=True)\n",
" f += -1j*delta/omega*bndc*grad(v)*n * ds(bnd,skeleton=True)\n",
" \n",
" with TaskManager():\n",
" a.Assemble()\n",
" f.Assemble()\n",
" return a,f"
]
},
{
Expand All @@ -256,16 +249,49 @@
"metadata": {},
"outputs": [],
"source": [
"mesh = Mesh(unit_square.GenerateMesh(maxh=.3))\n",
"order = 4\n",
"fes = L2(mesh, order=order, complex=True, dgjumps=True)\n",
"a,f = dghelmholtz(fes,omega,bndc=bndc)\n",
"PP = TrefftzHelmholtzEmb(fes)\n",
"PPT = PP.CreateTranspose()\n",
"with TaskManager():\n",
" TA = PPT@a.mat@PP\n",
" TU = TA.Inverse()*(PPT*f.vec)\n",
" tpgfu = GridFunction(fes)\n",
" tpgfu.vec.data = PP*TU\n",
"error = sqrt(Integrate((tpgfu-exact)*Conj(tpgfu-exact), mesh).real)\n",
" gfu = GridFunction(fes)\n",
" gfu.vec.data = PP*TU\n",
"error = sqrt(Integrate((gfu-exact)*Conj(gfu-exact), mesh).real)\n",
"print(\"error \",error)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "091833e6-f6d4-4ed6-8589-1ddbc47b01e7",
"metadata": {},
"outputs": [],
"source": [
"air = Circle((0.5, 0.5), 0.8).Face()\n",
"air.edges.name = 'outer'\n",
"scatterer = MoveTo(0.7, 0.3).Rectangle(0.05, 0.4).Face()\n",
"scatterer.edges.name = 'scat'\n",
"geo = OCCGeometry(air - scatterer, dim=2)\n",
"mesh = Mesh(geo.GenerateMesh(maxh=0.05))\n",
"omega = 25\n",
"rhs = 3e2*exp(-(40**2)*((x-0.5)*(x-0.5) + (y-0.5)*(y-0.5)))\n",
"\n",
"fes = L2(mesh, order=order, complex=True, dgjumps=True)\n",
"a,f = dghelmholtz(fes,omega,rhs=rhs,bnd=\"outer\")\n",
"PP = TrefftzHelmholtzEmb(fes,omega)\n",
"PPT = PP.CreateTranspose()\n",
"with TaskManager():\n",
" TA = PPT@a.mat@PP\n",
" TU = TA.Inverse()*(PPT*f.vec)\n",
" gfu = GridFunction(fes)\n",
" gfu.vec.data = PP*TU\n",
"\n",
"Draw(gfu, mesh, animate_complex=True,deformation=True);"
]
}
],
"metadata": {
Expand All @@ -284,7 +310,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.11.3"
"version": "3.12.3"
}
},
"nbformat": 4,
Expand Down
4 changes: 2 additions & 2 deletions docs/notebooks/helmholtz.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@
"\\begin{align*}\n",
" \\begin{cases}\n",
" -\\Delta u - \\omega^2 u= 0 &\\text{ in } \\dom, \\\\\n",
" \\frac{\\partial u}{\\partial \\nx} + i u = g &\\text{ on } \\partial \\dom.\n",
" \\frac{\\partial u}{\\partial \\nx} + i\\omega u = g &\\text{ on } \\partial \\dom.\n",
" \\end{cases}\n",
"\\end{align*}\n",
"$$\n"
Expand Down Expand Up @@ -173,7 +173,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.11.7"
"version": "3.12.3"
}
},
"nbformat": 4,
Expand Down
51 changes: 12 additions & 39 deletions docs/notebooks/qtelliptic.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@
"id": "1ffc608a-0a36-49c0-9a86-646e3b82ca20",
"metadata": {},
"source": [
"# Quasi-Treffz DG: diffusion-advection-reaction \n",
"# Quasi-Trefftz DG: diffusion-advection-reaction \n",
"Consider the following diffusion-advection-reaction BVP\n",
"$$\n",
"\\begin{cases}\n",
Expand All @@ -28,7 +28,8 @@
"source": [
"from ngsolve import *\n",
"from ngstrefftz import *\n",
"from netgen.occ import *"
"from netgen.occ import *\n",
"SetNumThreads(4)"
]
},
{
Expand Down Expand Up @@ -63,9 +64,7 @@
"h = 0.03; \n",
"order = 3;\n",
"mesh = Mesh(unit_square.GenerateMesh(maxh=h))\n",
"SetNumThreads(4)\n",
"with TaskManager():\n",
" fes = trefftzfespace(mesh,order=order,eq=\"qtelliptic\")"
"fes = trefftzfespace(mesh,order=order,eq=\"qtelliptic\")"
]
},
{
Expand Down Expand Up @@ -103,7 +102,7 @@
"id": "0e26cef3-2c9c-4924-b0b3-45d59eb3cfd4",
"metadata": {},
"source": [
"To set the coefficients for the construction of the quasi-Trefftz basis functions of $ \\mathbb{Q\\!T}^p_0(\\mathcal{T}_h)$ use "
"To set the coefficients for the construction of the quasi-Trefftz basis functions of $\\mathbb{Q\\!T}^p_0(\\mathcal{T}_h)$ use "
]
},
{
Expand Down Expand Up @@ -243,26 +242,8 @@
"gfu = GridFunction(fes)\n",
"\n",
"with TaskManager():\n",
" gfu.vec.data = a.mat.Inverse()*f.vec"
]
},
{
"cell_type": "markdown",
"id": "0cbc4795-cb36-4129-8d76-c96091d8c443",
"metadata": {},
"source": [
"We can calculate the $L^2(\\Omega)$-error:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ff7a8109-3657-40a5-bfe1-ed08b8efbd5d",
"metadata": {},
"outputs": [],
"source": [
"error = sqrt(Integrate((gfu-exact)**2, mesh))\n",
"print(\"quasi-Trefftz DG error:\",error)"
" gfu.vec.data = a.mat.Inverse()*f.vec\n",
"print(\"quasi-Trefftz DG error:\",sqrt(Integrate((gfu-exact)**2, mesh)))"
]
},
{
Expand Down Expand Up @@ -327,15 +308,7 @@
"gfu = GridFunction(fes)\n",
"with TaskManager():\n",
" gfu.vec.data = a.mat.Inverse() * f.vec\n",
"Draw(gfu,mesh,\"gfu\")"
]
},
{
"cell_type": "markdown",
"id": "128e02a2-5115-476c-8dbc-e83eeb8ba48d",
"metadata": {},
"source": [
"To see the 3D surface click `Open Controls` and select `deformation`."
"Draw(gfu,mesh,\"gfu\",deformation=CF((x,y,gfu)))"
]
},
{
Expand All @@ -354,6 +327,7 @@
"outputs": [],
"source": [
"h = 0.05\n",
"order = 3\n",
"wp = WorkPlane().RectangleC(2,2) \\\n",
" .Circle(0,0,0.25).Reverse() \n",
"geo = wp.Face()\n",
Expand Down Expand Up @@ -386,7 +360,6 @@
"rhs = 0\n",
"uf = None\n",
"\n",
"order = 3\n",
"with TaskManager():\n",
" fes = trefftzfespace(mesh,order=order,eq=\"qtelliptic\")\n",
" fes.SetCoeff(K,beta,sigma)\n",
Expand Down Expand Up @@ -430,7 +403,7 @@
"metadata": {},
"outputs": [],
"source": [
"nu = 5*10**(-3)\n",
"nu = 10**(-2)\n",
"K = CF((nu,0,0,nu),dims=(2,2))\n",
"beta = CF((-y,x))\n",
"sigma = 0\n",
Expand All @@ -449,7 +422,7 @@
"gfu = GridFunction(fes)\n",
"with TaskManager():\n",
" gfu.vec.data = a.mat.Inverse() * f.vec\n",
"Draw(gfu,mesh,\"gfu\")"
"Draw(gfu,mesh,\"gfu\",deformation=CF((x,y,gfu)))"
]
}
],
Expand All @@ -469,7 +442,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.12.2"
"version": "3.12.3"
}
},
"nbformat": 4,
Expand Down
2 changes: 1 addition & 1 deletion docs/notebooks/qtwave.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@
"id": "5e251ded",
"metadata": {},
"source": [
"# Quasi-Trefftz DG "
"# Quasi-Trefftz DG: Wave equation"
]
},
{
Expand Down
2 changes: 1 addition & 1 deletion docs/notebooks/twave.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -281,7 +281,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.11.3"
"version": "3.12.3"
}
},
"nbformat": 4,
Expand Down

0 comments on commit c0ab962

Please sign in to comment.