|
21 | 21 | "\\end{align*}\n", |
22 | 22 | "for the specific choice\n", |
23 | 23 | "\\begin{align*}\n", |
24 | | - " f(x, y) = 8\\pi^2\\sin(2\\pi x)\\sin(2\\pi y)\n", |
| 24 | + " f(x, y) = [ (-x^4 + x^3 + (-y^2 + \\pi^2)x^2 + (y^2 - 4y - \\pi^2)x + 2y - 2) \\sin(\\pi y) + (2\\pi x^2(1-x)) \\cos(\\pi y)] e^{xy}\n", |
25 | 25 | "\\end{align*}\n", |
26 | 26 | "\n", |
27 | 27 | "## The Variational Formulation\n", |
|
63 | 63 | "phi, psi = elements_of(V, names='phi, psi')\n", |
64 | 64 | "\n", |
65 | 65 | "# bilinear form\n", |
66 | | - "a = BilinearForm((phi, psi), integral(domain, dot(grad(phi), grad(psi))))\n", |
| 66 | + "a = BilinearForm((phi, psi), integral(domain, phi*psi))\n", |
67 | 67 | "\n", |
68 | 68 | "# linear form\n", |
69 | | - "from sympy import pi, sin\n", |
| 69 | + "from sympy import pi, sin, cos, exp\n", |
70 | 70 | "x,y = domain.coordinates\n", |
71 | | - "f = 8 * (pi**2) * sin(2 * pi * x) * sin(2 * pi * y)\n", |
| 71 | + "f = ( (-x**4 + x**3 + (-y**2 + pi**2)*x**2 + (y**2 - 4*y - pi**2)*x + 2*y - 2) * sin(pi*y) + (2*pi*x**2*(1-x)) * cos(pi*y) ) * exp(x*y)\n", |
72 | 72 | "\n", |
73 | 73 | "l = LinearForm(psi, integral(domain, f*psi))" |
74 | 74 | ] |
|
167 | 167 | "metadata": {}, |
168 | 168 | "outputs": [], |
169 | 169 | "source": [ |
170 | | - "from psydac.linalg.solvers import inverse\n", |
| 170 | + "import time\n", |
| 171 | + "\n", |
| 172 | + "from psydac.linalg.solvers import inverse\n", |
171 | 173 | "\n", |
172 | 174 | "tol = 1e-12\n", |
173 | 175 | "maxiter = 1000\n", |
174 | 176 | "\n", |
175 | | - "phi_h = inverse(A_bc, 'cg', tol=tol, maxiter=maxiter) @ f_bc" |
| 177 | + "A_bc_inv = inverse(A_bc, 'cg', tol=tol, maxiter=maxiter)\n", |
| 178 | + "\n", |
| 179 | + "t0 = time.time()\n", |
| 180 | + "phi_h = A_bc_inv @ f_bc\n", |
| 181 | + "t1 = time.time()" |
176 | 182 | ] |
177 | 183 | }, |
178 | 184 | { |
|
186 | 192 | "In this example, the analytical solution is given by\n", |
187 | 193 | "\n", |
188 | 194 | "$$\n", |
189 | | - "phi_{ex}(x, y) = \\sin(2\\pi x)\\sin(2\\pi y)\n", |
| 195 | + "phi_{ex}(x, y) = x(x-1)\\sin(\\pi y)e^{xy}\n", |
190 | 196 | "$$" |
191 | 197 | ] |
192 | 198 | }, |
|
207 | 213 | "\n", |
208 | 214 | "from sympde.expr import Norm, SemiNorm\n", |
209 | 215 | "\n", |
210 | | - "phi_ex = sin(2*pi*x) * sin(2*pi*y)\n", |
| 216 | + "phi_ex = x*(x-1)*sin(pi*y)*exp(x*y)\n", |
211 | 217 | "phi_h_FemField = FemField(V_h, phi_h)\n", |
212 | 218 | "\n", |
213 | 219 | "error = phi_ex - phi\n", |
|
219 | 225 | "l2norm_h = discretize(l2norm, domain_h, V_h, backend=backend)\n", |
220 | 226 | "\n", |
221 | 227 | "# assemble the norm\n", |
222 | | - "l2_error = l2norm_h.assemble(phi=phi_h_FemField)\n", |
223 | | - "\n", |
224 | | - "# print the result\n", |
225 | | - "print(l2_error)" |
| 228 | + "l2_error = l2norm_h.assemble(phi=phi_h_FemField)" |
226 | 229 | ] |
227 | 230 | }, |
228 | 231 | { |
|
245 | 248 | "h1norm_h = discretize(h1norm, domain_h, V_h)\n", |
246 | 249 | "\n", |
247 | 250 | "# assemble the norm\n", |
248 | | - "h1_error = h1norm_h.assemble(phi=phi_h_FemField)\n", |
249 | | - "\n", |
250 | | - "# print the result\n", |
251 | | - "print(h1_error)" |
| 251 | + "h1semi_error = h1norm_h.assemble(phi=phi_h_FemField)" |
252 | 252 | ] |
253 | 253 | }, |
254 | 254 | { |
|
274 | 274 | "h1_error = h1norm_h.assemble(phi=phi_h_FemField)\n", |
275 | 275 | "\n", |
276 | 276 | "# print the result\n", |
277 | | - "print(h1_error)" |
| 277 | + "print( '> Grid :: [{ne1},{ne2}]'.format( ne1=ncells[0], ne2=ncells[1]) )\n", |
| 278 | + "print( '> Degree :: [{p1},{p2}]' .format( p1=degree[0], p2=degree[1] ) )\n", |
| 279 | + "print( '> CG info :: ',A_bc_inv.get_info() )\n", |
| 280 | + "print( '> L2 error :: {:.2e}'.format( l2_error ) )\n", |
| 281 | + "print( '> H1-Semi error :: {:.2e}'.format( h1semi_error ) )\n", |
| 282 | + "print( '> H1 error :: {:.2e}'.format( h1_error ) )\n", |
| 283 | + "print( '' )\n", |
| 284 | + "print( '> Solution time :: {:.3g}'.format( t1-t0 ) )" |
278 | 285 | ] |
279 | 286 | }, |
280 | 287 | { |
|
0 commit comments