Skip to content

Commit 6dcdb41

Browse files
jowezarekjowezarek
committed
add sheet2 exercise 2.3
1 parent cada87e commit 6dcdb41

File tree

2 files changed

+308
-8
lines changed

2 files changed

+308
-8
lines changed

exercises/time_harmonic_maxwell_2d.ipynb renamed to exercises/time_harmonic_maxwell_2d_I.ipynb

Lines changed: 9 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -12,6 +12,7 @@
1212
"metadata": {},
1313
"source": [
1414
"# 2D Solver for the time-harmonic Maxwell equation on the Unit Square\n",
15+
"## ... using Psydac's bilinear form interface\n",
1516
"\n",
1617
"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",
1718
"\n",
@@ -100,7 +101,7 @@
100101
"from psydac.api.discretization import discretize\n",
101102
"from psydac.api.settings import PSYDAC_BACKEND_GPYCCEL\n",
102103
"\n",
103-
"backend = PSYDAC_BACKEND_GPYCCEL"
104+
"backend = PSYDAC_BACKEND_GPYCCEL"
104105
]
105106
},
106107
{
@@ -180,7 +181,7 @@
180181
"\n",
181182
"from psydac.linalg.solvers import inverse\n",
182183
"\n",
183-
"tol = 1e-12\n",
184+
"tol = 1e-9\n",
184185
"maxiter = 1000\n",
185186
"\n",
186187
"A_bc_inv = inverse(A_bc, 'cg', tol=tol, maxiter=maxiter)\n",
@@ -217,19 +218,19 @@
217218
"\n",
218219
"from sympde.expr import Norm\n",
219220
"\n",
220-
"E_ex = Tuple(sin(pi*y), sin(pi*x)*cos(pi*y))\n",
221-
"E_h_FemField = FemField(V_h, E_h)\n",
221+
"E_ex = Tuple(sin(pi*y), sin(pi*x)*cos(pi*y))\n",
222+
"E_h_FemField = FemField(V_h, E_h)\n",
222223
"\n",
223-
"error = Matrix([E[0] - E_ex[0], E[1] - E_ex[1]])\n",
224+
"error = Matrix([E[0] - E_ex[0], E[1] - E_ex[1]])\n",
224225
"\n",
225226
"# create the formal Norm object\n",
226-
"l2norm = Norm(error, domain, kind='l2')\n",
227+
"l2norm = Norm(error, domain, kind='l2')\n",
227228
"\n",
228229
"# discretize the norm\n",
229-
"l2norm_h = discretize(l2norm, domain_h, V_h, backend=backend)\n",
230+
"l2norm_h = discretize(l2norm, domain_h, V_h, backend=backend)\n",
230231
"\n",
231232
"# assemble the norm\n",
232-
"l2_error = l2norm_h.assemble(E=E_h_FemField)\n",
233+
"l2_error = l2norm_h.assemble(E=E_h_FemField)\n",
233234
"\n",
234235
"# print the result\n",
235236
"print( '> Grid :: [{ne1},{ne2}]'.format( ne1=ncells[0], ne2=ncells[1]) )\n",

‎exercises/time_harmonic_maxwell_2d_II.ipynb‎

Lines changed: 299 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,299 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "markdown",
5+
"metadata": {},
6+
"source": [
7+
"**Note: This notebook must be changed by the student. As it is now, it does not approximate the solution to the Poisson problem!**"
8+
]
9+
},
10+
{
11+
"cell_type": "markdown",
12+
"metadata": {},
13+
"source": [
14+
"# 2D Solver for the time-harmonic Maxwell equation on the Unit Square\n",
15+
"## ... using Psydac's de Rham interface\n",
16+
"\n",
17+
"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",
18+
"\n",
19+
"\\begin{align*}\n",
20+
" \\boldsymbol{curl}curl\\boldsymbol{E} - \\omega^2\\boldsymbol{E} &= \\boldsymbol{F} \\quad \\text{ in }\\Omega=]0,1[^2, \\\\\n",
21+
" \\boldsymbol{n}\\times\\boldsymbol{E} &= 0 \\quad \\text{ on }\\partial\\Omega,\n",
22+
"\\end{align*}\n",
23+
"for the specific choice\n",
24+
"\\begin{align*}\n",
25+
" \\omega &= 1.5 \\\\\n",
26+
" \\alpha &= -\\omega^2 \\\\\n",
27+
" \\boldsymbol{F}(x, y) &= \\left(\\begin{matrix}\n",
28+
" (\\alpha+\\pi^2)\\sin(\\pi y) - \\pi^2\\cos(\\pi x)\\sin(\\pi y) \\\\\n",
29+
" (\\alpha+\\pi^2)\\sin(\\pi x)\\cos(\\pi y)\n",
30+
" \\end{matrix}\\right)\n",
31+
"\\end{align*}\n",
32+
"\n",
33+
"## The Variational Formulation\n",
34+
"\n",
35+
"The corresponding variational formulation reads\n",
36+
"\n",
37+
"\\begin{align*}\n",
38+
" \\text{Find }\\boldsymbol{E}\\in V\\coloneqq H_0(curl;\\Omega)\\text{ such that } \\\\\n",
39+
" a(\\boldsymbol{E}, \\boldsymbol{v}) = l(\\boldsymbol{v})\\quad\\forall\\ \\boldsymbol{v}\\in V\n",
40+
"\\end{align*}\n",
41+
"\n",
42+
"where \n",
43+
"- $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",
44+
"- $l(\\boldsymbol{v}) := \\int_{\\Omega} \\boldsymbol{F}\\cdot\\boldsymbol{v} ~ d\\Omega$."
45+
]
46+
},
47+
{
48+
"cell_type": "markdown",
49+
"metadata": {},
50+
"source": [
51+
"## Discrete Model using de Rham objects"
52+
]
53+
},
54+
{
55+
"cell_type": "code",
56+
"execution_count": null,
57+
"metadata": {},
58+
"outputs": [],
59+
"source": [
60+
"from sympde.calculus import dot\n",
61+
"from sympde.topology import elements_of, Square, Derham\n",
62+
"from sympde.expr.expr import LinearForm, BilinearForm, integral\n",
63+
"\n",
64+
"domain = Square('S', bounds1=(0, 1), bounds2=(0, 1))\n",
65+
"derham = Derham(domain, sequence=['h1', 'hcurl', 'l2'])\n",
66+
"\n",
67+
"V1 = derham.V1\n",
68+
"V2 = derham.V2\n",
69+
"\n",
70+
"u1, v1 = elements_of(V1, names='u1, v1')\n",
71+
"u2, v2 = elements_of(V2, names='u2, v2')\n",
72+
"\n",
73+
"# bilinear forms corresponding to V1 and V2 mass matrices\n",
74+
"m1 = BilinearForm((u1, v1), integral(domain, dot(u1, v1)))\n",
75+
"m2 = BilinearForm((u2, v2), integral(domain, u2 * v2))\n",
76+
"\n",
77+
"# linear form corresponding to the rhs\n",
78+
"from sympy import pi, sin, cos, Tuple, Matrix\n",
79+
"omega = 1.5\n",
80+
"alpha = -omega**2\n",
81+
"x,y = domain.coordinates\n",
82+
"F = Tuple( (alpha + pi**2) * sin(pi*y) - pi**2 * sin(pi*y) * cos(pi*x),\n",
83+
" (alpha + pi**2) * sin(pi*x) * cos(pi*y) )\n",
84+
"\n",
85+
"expr = dot(F, v1)\n",
86+
"l = LinearForm(v1, integral(domain, expr))"
87+
]
88+
},
89+
{
90+
"cell_type": "code",
91+
"execution_count": null,
92+
"metadata": {},
93+
"outputs": [],
94+
"source": [
95+
"from psydac.api.discretization import discretize\n",
96+
"from psydac.api.settings import PSYDAC_BACKEND_GPYCCEL\n",
97+
"\n",
98+
"backend = PSYDAC_BACKEND_GPYCCEL"
99+
]
100+
},
101+
{
102+
"cell_type": "code",
103+
"execution_count": null,
104+
"metadata": {},
105+
"outputs": [],
106+
"source": [
107+
"ncells = [16, 16] # Bspline cells\n",
108+
"degree = [3, 3] # Bspline degree"
109+
]
110+
},
111+
{
112+
"cell_type": "code",
113+
"execution_count": null,
114+
"metadata": {},
115+
"outputs": [],
116+
"source": [
117+
"domain_h = discretize(domain, ncells=ncells, periodic=[False, False])\n",
118+
"derham_h = discretize(derham, domain_h, degree=degree)\n",
119+
"\n",
120+
"V1_h = derham_h.V1\n",
121+
"V2_h = derham_h.V2\n",
122+
"\n",
123+
"# Exterior Derivative operators (grad and curl)\n",
124+
"G, C = derham_h.derivatives_as_matrices\n",
125+
"\n",
126+
"# Mass matrices\n",
127+
"m1_h = discretize(m1, domain_h, (V1_h, V1_h), backend=backend)\n",
128+
"m2_h = discretize(m2, domain_h, (V2_h, V2_h), backend=backend)\n",
129+
"\n",
130+
"M1 = m1_h.assemble()\n",
131+
"M2 = m2_h.assemble()\n",
132+
"\n",
133+
"# System Matrix A\n",
134+
"A = M1\n",
135+
"\n",
136+
"# rhs vector f\n",
137+
"l_h = discretize(l, domain_h, V1_h, backend=backend)\n",
138+
"f = l_h.assemble()"
139+
]
140+
},
141+
{
142+
"cell_type": "markdown",
143+
"metadata": {},
144+
"source": [
145+
"## Boundary Conditions\n",
146+
"\n",
147+
"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$."
148+
]
149+
},
150+
{
151+
"cell_type": "code",
152+
"execution_count": null,
153+
"metadata": {},
154+
"outputs": [],
155+
"source": [
156+
"from utils import HcurlBoundaryProjector2D\n",
157+
"from psydac.linalg.basic import IdentityOperator\n",
158+
"\n",
159+
"P_0 = HcurlBoundaryProjector2D(V1, V1_h.vector_space)\n",
160+
"\n",
161+
"I_0 = IdentityOperator(V1_h.vector_space)\n",
162+
"P_Gamma = I_0 - P_0"
163+
]
164+
},
165+
{
166+
"cell_type": "code",
167+
"execution_count": null,
168+
"metadata": {},
169+
"outputs": [],
170+
"source": [
171+
"A_bc = P_0 @ A @ P_0 + P_Gamma\n",
172+
"f_bc = P_0 @ f"
173+
]
174+
},
175+
{
176+
"cell_type": "markdown",
177+
"metadata": {},
178+
"source": [
179+
"## Solving the PDE"
180+
]
181+
},
182+
{
183+
"cell_type": "code",
184+
"execution_count": null,
185+
"metadata": {},
186+
"outputs": [],
187+
"source": [
188+
"import time\n",
189+
"\n",
190+
"from psydac.linalg.solvers import inverse\n",
191+
"\n",
192+
"tol = 1e-9\n",
193+
"maxiter = 1000\n",
194+
"\n",
195+
"A_bc_inv = inverse(A_bc, 'cg', tol=tol, maxiter=maxiter)\n",
196+
"\n",
197+
"t0 = time.time()\n",
198+
"E_h = A_bc_inv @ f_bc\n",
199+
"t1 = time.time()"
200+
]
201+
},
202+
{
203+
"cell_type": "markdown",
204+
"metadata": {},
205+
"source": [
206+
"## Computing the error norm\n",
207+
"\n",
208+
"As the analytical solution is available, we want to compute the $L^2$ norm of the error.\n",
209+
"In this example, the analytical solution is given by\n",
210+
"\n",
211+
"$$\n",
212+
"\\boldsymbol{E}_{ex}(x, y) = \\left(\\begin{matrix} \n",
213+
" \\sin(\\pi y) \\\\\n",
214+
" \\sin(\\pi x)\\cos(\\pi y)\n",
215+
" \\end{matrix}\\right)\n",
216+
"$$"
217+
]
218+
},
219+
{
220+
"cell_type": "code",
221+
"execution_count": null,
222+
"metadata": {},
223+
"outputs": [],
224+
"source": [
225+
"from psydac.fem.basic import FemField\n",
226+
"\n",
227+
"from sympde.expr import Norm\n",
228+
"\n",
229+
"E_ex = Tuple(sin(pi*y), sin(pi*x)*cos(pi*y))\n",
230+
"E_h_FemField = FemField(V1_h, E_h)\n",
231+
"\n",
232+
"error = Matrix([u1[0] - E_ex[0], u1[1] - E_ex[1]])\n",
233+
"\n",
234+
"# create the formal Norm object\n",
235+
"l2norm = Norm(error, domain, kind='l2')\n",
236+
"\n",
237+
"# discretize the norm\n",
238+
"l2norm_h = discretize(l2norm, domain_h, V1_h, backend=backend)\n",
239+
"\n",
240+
"# assemble the norm\n",
241+
"l2_error = l2norm_h.assemble(u1=E_h_FemField)\n",
242+
"\n",
243+
"# print the result\n",
244+
"print( '> Grid :: [{ne1},{ne2}]'.format( ne1=ncells[0], ne2=ncells[1]) )\n",
245+
"print( '> Degree :: [{p1},{p2}]' .format( p1=degree[0], p2=degree[1] ) )\n",
246+
"print( '> CG info :: ',A_bc_inv.get_info() )\n",
247+
"print( '> L2 error :: {:.2e}'.format( l2_error ) )\n",
248+
"print( '' )\n",
249+
"print( '> Solution time :: {:.3g}'.format( t1-t0 ) )"
250+
]
251+
},
252+
{
253+
"cell_type": "markdown",
254+
"metadata": {},
255+
"source": [
256+
"## Visualization\n",
257+
"\n",
258+
"We plot the true solution $\\boldsymbol{E}_{ex}$, the approximate solution $\\boldsymbol{E}_h$ and the error function $|\\boldsymbol{E}_{ex} - \\boldsymbol{E}_h|$."
259+
]
260+
},
261+
{
262+
"cell_type": "code",
263+
"execution_count": null,
264+
"metadata": {},
265+
"outputs": [],
266+
"source": [
267+
"from sympy import lambdify\n",
268+
"\n",
269+
"from utils import plot\n",
270+
"\n",
271+
"E_ex_x = lambdify((x, y), E_ex[0])\n",
272+
"E_ex_y = lambdify((x, y), E_ex[1])\n",
273+
"E_h_x = E_h_FemField[0]\n",
274+
"E_h_y = E_h_FemField[1]\n",
275+
"error_x = lambda x, y: abs(E_ex_x(x, y) - E_h_x(x, y))\n",
276+
"error_y = lambda x, y: abs(E_ex_y(x, y) - E_h_y(x, y))\n",
277+
"\n",
278+
"plot(gridsize_x = 100, \n",
279+
" gridsize_y = 100, \n",
280+
" title = r'Approximation of Solution $\\boldsymbol{E}$, first component', \n",
281+
" funs = [E_ex_x, E_h_x, error_x], \n",
282+
" 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",
283+
" surface_plot = True\n",
284+
")\n",
285+
"\n",
286+
"plot(gridsize_x = 100,\n",
287+
" gridsize_y = 100,\n",
288+
" title = r'Approximation of Solution $\\boldsymbol{E}$, second component',\n",
289+
" funs = [E_ex_y, E_h_y, error_y],\n",
290+
" 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",
291+
" surface_plot = True\n",
292+
")"
293+
]
294+
}
295+
],
296+
"metadata": {},
297+
"nbformat": 4,
298+
"nbformat_minor": 2
299+
}

0 commit comments

Comments
 (0)