Skip to content

Commit 763405b

Browse files
jowezarekjowezarek
committed
add sheet1 exercise 1.3 iii)
1 parent f9537c0 commit 763405b

File tree

1 file changed

+304
-0
lines changed

1 file changed

+304
-0
lines changed

‎exercises/poisson_2d_I_naturalBC.ipynb‎

Lines changed: 304 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,304 @@
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 Poisson Solver on the Unit Square\n",
15+
"## ... using Psydac's bilinear form interface\n",
16+
"\n",
17+
"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",
18+
"\n",
19+
"\\begin{align*}\n",
20+
" -\\Delta\\phi &= f \\quad \\text{ in }\\Omega=]0,1[^2, \\\\\n",
21+
" \\frac{\\partial\\phi}{\\partial \\boldsymbol{n}} &= 0 \\quad \\text{ on }\\partial\\Omega,\n",
22+
"\\end{align*}\n",
23+
"for the specific choice\n",
24+
"\\begin{align*}\n",
25+
" f(x, y) = 8\\pi^2\\cos(2\\pi x)\\cos(2\\pi y)\n",
26+
"\\end{align*}\n",
27+
"\n",
28+
"## The Variational Formulation\n",
29+
"\n",
30+
"The corresponding variational formulation reads\n",
31+
"\n",
32+
"\\begin{align*}\n",
33+
" \\text{Find }\\phi\\in V\\coloneqq H^1(\\Omega)\\text{ such that } \\\\\n",
34+
" a(\\phi, \\psi) = l(\\psi)\\quad\\forall\\ \\psi\\in V\n",
35+
"\\end{align*}\n",
36+
"\n",
37+
"where \n",
38+
"- $a(\\phi,\\psi) \\coloneqq \\int_{\\Omega} \\nabla\\phi\\cdot\\nabla\\psi ~ d\\Omega$,\n",
39+
"- $l(\\psi) := \\int_{\\Omega} f\\psi ~ d\\Omega$."
40+
]
41+
},
42+
{
43+
"cell_type": "markdown",
44+
"metadata": {},
45+
"source": [
46+
"## Formal Model"
47+
]
48+
},
49+
{
50+
"cell_type": "code",
51+
"execution_count": null,
52+
"metadata": {},
53+
"outputs": [],
54+
"source": [
55+
"from sympde.calculus import dot, grad\n",
56+
"from sympde.expr import BilinearForm, LinearForm, integral\n",
57+
"from sympde.topology import elements_of, Square, ScalarFunctionSpace\n",
58+
"\n",
59+
"domain = Square('S', bounds1=(0, 1), bounds2=(0, 1))\n",
60+
"\n",
61+
"V = ScalarFunctionSpace('V', domain, kind='h1')\n",
62+
"\n",
63+
"phi, psi = elements_of(V, names='phi, psi')\n",
64+
"\n",
65+
"# bilinear form\n",
66+
"a = BilinearForm((phi, psi), integral(domain, dot(grad(phi), grad(psi))))\n",
67+
"\n",
68+
"# linear form\n",
69+
"from sympy import pi, sin, cos, exp\n",
70+
"x,y = domain.coordinates\n",
71+
"f = 8 * pi**2 * cos(2*pi*x) * cos(2*pi*y)\n",
72+
"\n",
73+
"l = LinearForm(psi, integral(domain, f*psi))"
74+
]
75+
},
76+
{
77+
"cell_type": "markdown",
78+
"metadata": {},
79+
"source": [
80+
"## Discretization\n",
81+
"\n",
82+
"We will use Psydac to discretize the problem."
83+
]
84+
},
85+
{
86+
"cell_type": "code",
87+
"execution_count": null,
88+
"metadata": {},
89+
"outputs": [],
90+
"source": [
91+
"from psydac.api.discretization import discretize\n",
92+
"from psydac.api.settings import PSYDAC_BACKEND_GPYCCEL\n",
93+
"\n",
94+
"backend = PSYDAC_BACKEND_GPYCCEL"
95+
]
96+
},
97+
{
98+
"cell_type": "code",
99+
"execution_count": null,
100+
"metadata": {},
101+
"outputs": [],
102+
"source": [
103+
"ncells = [16, 16] # Bspline cells\n",
104+
"degree = [3, 3] # Bspline degree"
105+
]
106+
},
107+
{
108+
"cell_type": "code",
109+
"execution_count": null,
110+
"metadata": {},
111+
"outputs": [],
112+
"source": [
113+
"domain_h = discretize(domain, ncells=ncells, periodic=[False, False])\n",
114+
"V_h = discretize(V, domain_h, degree=degree)\n",
115+
"\n",
116+
"a_h = discretize(a, domain_h, (V_h, V_h), backend=backend)\n",
117+
"l_h = discretize(l, domain_h, V_h, backend=backend)"
118+
]
119+
},
120+
{
121+
"cell_type": "markdown",
122+
"metadata": {},
123+
"source": [
124+
"## Boundary Conditions\n",
125+
"\n",
126+
"We must not take care of boundary conditions. They are included **naturally** in the variational formulation."
127+
]
128+
},
129+
{
130+
"cell_type": "code",
131+
"execution_count": null,
132+
"metadata": {},
133+
"outputs": [],
134+
"source": [
135+
"A = a_h.assemble()\n",
136+
"f = l_h.assemble()"
137+
]
138+
},
139+
{
140+
"cell_type": "markdown",
141+
"metadata": {},
142+
"source": [
143+
"## Solving the PDE"
144+
]
145+
},
146+
{
147+
"cell_type": "code",
148+
"execution_count": null,
149+
"metadata": {},
150+
"outputs": [],
151+
"source": [
152+
"import time\n",
153+
"\n",
154+
"from psydac.linalg.solvers import inverse\n",
155+
"\n",
156+
"tol = 1e-12\n",
157+
"maxiter = 1000\n",
158+
"\n",
159+
"A_bc_inv = inverse(A, 'cg', tol=tol, maxiter=maxiter)\n",
160+
"\n",
161+
"t0 = time.time()\n",
162+
"phi_h = A_bc_inv @ f\n",
163+
"t1 = time.time()"
164+
]
165+
},
166+
{
167+
"cell_type": "markdown",
168+
"metadata": {},
169+
"source": [
170+
"## Computing the error norm\n",
171+
"\n",
172+
"When the analytical solution is available, you might be interested in computing the $L^2$ norm or $H^1_0$ semi-norm.\n",
173+
"SymPDE allows you to do so, by creating the **Norm** object.\n",
174+
"In this example, the analytical solution is given by\n",
175+
"\n",
176+
"$$\n",
177+
"phi_{ex}(x, y) = \\cos(2\\pi x)\\cos(2\\pi y)\n",
178+
"$$"
179+
]
180+
},
181+
{
182+
"cell_type": "markdown",
183+
"metadata": {},
184+
"source": [
185+
"### Computing the $L^2$ norm"
186+
]
187+
},
188+
{
189+
"cell_type": "code",
190+
"execution_count": null,
191+
"metadata": {},
192+
"outputs": [],
193+
"source": [
194+
"from psydac.fem.basic import FemField\n",
195+
"\n",
196+
"from sympde.expr import Norm, SemiNorm\n",
197+
"\n",
198+
"phi_ex = cos(2*pi*x) * cos(2*pi*y)\n",
199+
"phi_h_FemField = FemField(V_h, phi_h)\n",
200+
"\n",
201+
"error = phi_ex - phi\n",
202+
"\n",
203+
"# create the formal Norm object\n",
204+
"l2norm = Norm(error, domain, kind='l2')\n",
205+
"\n",
206+
"# discretize the norm\n",
207+
"l2norm_h = discretize(l2norm, domain_h, V_h, backend=backend)\n",
208+
"\n",
209+
"# assemble the norm\n",
210+
"l2_error = l2norm_h.assemble(phi=phi_h_FemField)"
211+
]
212+
},
213+
{
214+
"cell_type": "markdown",
215+
"metadata": {},
216+
"source": [
217+
"### Computing the $H^1$ semi-norm"
218+
]
219+
},
220+
{
221+
"cell_type": "code",
222+
"execution_count": null,
223+
"metadata": {},
224+
"outputs": [],
225+
"source": [
226+
"# create the formal Norm object\n",
227+
"h1norm = SemiNorm(error, domain, kind='h1')\n",
228+
"\n",
229+
"# discretize the norm\n",
230+
"h1norm_h = discretize(h1norm, domain_h, V_h)\n",
231+
"\n",
232+
"# assemble the norm\n",
233+
"h1semi_error = h1norm_h.assemble(phi=phi_h_FemField)"
234+
]
235+
},
236+
{
237+
"cell_type": "markdown",
238+
"metadata": {},
239+
"source": [
240+
"### Computing the $H^1$ norm"
241+
]
242+
},
243+
{
244+
"cell_type": "code",
245+
"execution_count": null,
246+
"metadata": {},
247+
"outputs": [],
248+
"source": [
249+
"# create the formal Norm object\n",
250+
"h1norm = Norm(error, domain, kind='h1')\n",
251+
"\n",
252+
"# discretize the norm\n",
253+
"h1norm_h = discretize(h1norm, domain_h, V_h)\n",
254+
"\n",
255+
"# assemble the norm\n",
256+
"h1_error = h1norm_h.assemble(phi=phi_h_FemField)\n",
257+
"\n",
258+
"# print the result\n",
259+
"print( '> Grid :: [{ne1},{ne2}]'.format( ne1=ncells[0], ne2=ncells[1]) )\n",
260+
"print( '> Degree :: [{p1},{p2}]' .format( p1=degree[0], p2=degree[1] ) )\n",
261+
"print( '> CG info :: ',A_bc_inv.get_info() )\n",
262+
"print( '> L2 error :: {:.2e}'.format( l2_error ) )\n",
263+
"print( '> H1-Semi error :: {:.2e}'.format( h1semi_error ) )\n",
264+
"print( '> H1 error :: {:.2e}'.format( h1_error ) )\n",
265+
"print( '' )\n",
266+
"print( '> Solution time :: {:.3g}'.format( t1-t0 ) )"
267+
]
268+
},
269+
{
270+
"cell_type": "markdown",
271+
"metadata": {},
272+
"source": [
273+
"## Visualization\n",
274+
"\n",
275+
"We plot the true solution $\\phi_{ex}$, the approximate solution $\\phi_h$ and the error function $|\\phi_{ex} - \\phi_h|$."
276+
]
277+
},
278+
{
279+
"cell_type": "code",
280+
"execution_count": null,
281+
"metadata": {},
282+
"outputs": [],
283+
"source": [
284+
"from sympy import lambdify\n",
285+
"\n",
286+
"from utils import plot\n",
287+
"\n",
288+
"phi_ex_fun = lambdify(domain.coordinates, phi_ex)\n",
289+
"error = lambda x, y: abs(phi_ex_fun(x, y) - phi_h_FemField(x, y))\n",
290+
"\n",
291+
"plot(gridsize_x = 100, \n",
292+
" gridsize_y = 100, \n",
293+
" title = r'Approximation of Solution $\\phi$', \n",
294+
" funs = [phi_ex_fun, phi_h_FemField, error], \n",
295+
" titles = [r'$\\phi_{ex}(x,y)$', r'$\\phi_h(x,y)$', r'$|(\\phi_{ex}-\\phi_h)(x,y)|$'],\n",
296+
" surface_plot = True\n",
297+
")"
298+
]
299+
}
300+
],
301+
"metadata": {},
302+
"nbformat": 4,
303+
"nbformat_minor": 4
304+
}

0 commit comments

Comments
 (0)