Coverage for pySDC/implementations/problem_classes/AllenCahn_2D_FD.py: 96%
201 statements
« prev ^ index » next coverage.py v7.6.9, created at 2024-12-20 14:51 +0000
1import numpy as np
2import scipy.sparse as sp
3from scipy.sparse.linalg import cg
5from pySDC.core.errors import ParameterError, ProblemError
6from pySDC.core.problem import Problem, WorkCounter
7from pySDC.helpers import problem_helper
8from pySDC.implementations.datatype_classes.mesh import mesh, imex_mesh, comp2_mesh
11# http://www.personal.psu.edu/qud2/Res/Pre/dz09sisc.pdf
14# noinspection PyUnusedLocal
15class allencahn_fullyimplicit(Problem):
16 r"""
17 Example implementing the two-dimensional Allen-Cahn equation with periodic boundary conditions :math:`u \in [-1, 1]^2`
19 .. math::
20 \frac{\partial u}{\partial t} = \Delta u + \frac{1}{\varepsilon^2} u (1 - u^\nu)
22 for constant parameter :math:`\nu`. Initial condition are circles of the form
24 .. math::
25 u({\bf x}, 0) = \tanh\left(\frac{r - \sqrt{x_i^2 + y_j^2}}{\sqrt{2}\varepsilon}\right)
27 for :math:`i, j=0,..,N-1`, where :math:`N` is the number of spatial grid points. For time-stepping, the problem is
28 treated *fully-implicitly*, i.e., the nonlinear system is solved by Newton.
30 Parameters
31 ----------
32 nvars : tuple of int, optional
33 Number of unknowns in the problem, e.g. ``nvars=(128, 128)``.
34 nu : float, optional
35 Problem parameter :math:`\nu`.
36 eps : float, optional
37 Scaling parameter :math:`\varepsilon`.
38 newton_maxiter : int, optional
39 Maximum number of iterations for the Newton solver.
40 newton_tol : float, optional
41 Tolerance for Newton's method to terminate.
42 lin_tol : float, optional
43 Tolerance for linear solver to terminate.
44 lin_maxiter : int, optional
45 Maximum number of iterations for the linear solver.
46 radius : float, optional
47 Radius of the circles.
48 order : int, optional
49 Order of the finite difference matrix.
51 Attributes
52 ----------
53 A : scipy.spdiags
54 Second-order FD discretization of the 2D laplace operator.
55 dx : float
56 Distance between two spatial nodes (same for both directions).
57 xvalues : np.1darray
58 Spatial grid points, here both dimensions have the same grid points.
59 newton_itercount : int
60 Number of iterations of Newton solver.
61 lin_itercount
62 Number of iterations of linear solver.
63 newton_ncalls : int
64 Number of calls of Newton solver.
65 lin_ncalls : int
66 Number of calls of linear solver.
67 """
69 dtype_u = mesh
70 dtype_f = mesh
72 def __init__(
73 self,
74 nvars=(128, 128),
75 nu=2,
76 eps=0.04,
77 newton_maxiter=200,
78 newton_tol=1e-12,
79 lin_tol=1e-8,
80 lin_maxiter=100,
81 inexact_linear_ratio=None,
82 radius=0.25,
83 order=2,
84 ):
85 """Initialization routine"""
86 # we assert that nvars looks very particular here.. this will be necessary for coarsening in space later on
87 if len(nvars) != 2:
88 raise ProblemError('this is a 2d example, got %s' % nvars)
89 if nvars[0] != nvars[1]:
90 raise ProblemError('need a square domain, got %s' % nvars)
91 if nvars[0] % 2 != 0:
92 raise ProblemError('the setup requires nvars = 2^p per dimension')
94 # invoke super init, passing number of dofs, dtype_u and dtype_f
95 super().__init__((nvars, None, np.dtype('float64')))
96 self._makeAttributeAndRegister(
97 'nvars',
98 'nu',
99 'eps',
100 'radius',
101 'order',
102 localVars=locals(),
103 readOnly=True,
104 )
105 self._makeAttributeAndRegister(
106 'newton_maxiter',
107 'newton_tol',
108 'lin_tol',
109 'lin_maxiter',
110 'inexact_linear_ratio',
111 localVars=locals(),
112 readOnly=False,
113 )
115 # compute dx and get discretization matrix A
116 self.dx = 1.0 / self.nvars[0]
117 self.A, _ = problem_helper.get_finite_difference_matrix(
118 derivative=2,
119 order=self.order,
120 stencil_type='center',
121 dx=self.dx,
122 size=self.nvars[0],
123 dim=2,
124 bc='periodic',
125 )
126 self.xvalues = np.array([i * self.dx - 0.5 for i in range(self.nvars[0])])
128 self.newton_itercount = 0
129 self.lin_itercount = 0
130 self.newton_ncalls = 0
131 self.lin_ncalls = 0
133 self.work_counters['newton'] = WorkCounter()
134 self.work_counters['rhs'] = WorkCounter()
135 self.work_counters['linear'] = WorkCounter()
137 # noinspection PyTypeChecker
138 def solve_system(self, rhs, factor, u0, t):
139 """
140 Simple Newton solver.
142 Parameters
143 ----------
144 rhs : dtype_f
145 Right-hand side for the nonlinear system
146 factor : float
147 Abbrev. for the node-to-node stepsize (or any other factor required).
148 u0 : dtype_u
149 Initial guess for the iterative solver.
150 t : float
151 Current time (required here for the BC).
153 Returns
154 -------
155 me : dtype_u
156 The solution as mesh.
157 """
159 u = self.dtype_u(u0).flatten()
160 z = self.dtype_u(self.init, val=0.0).flatten()
161 nu = self.nu
162 eps2 = self.eps**2
164 Id = sp.eye(self.nvars[0] * self.nvars[1])
166 # start newton iteration
167 n = 0
168 res = 99
169 while n < self.newton_maxiter:
170 # form the function g with g(u) = 0
171 g = u - factor * (self.A.dot(u) + 1.0 / eps2 * u * (1.0 - u**nu)) - rhs.flatten()
173 # if g is close to 0, then we are done
174 res = np.linalg.norm(g, np.inf)
176 # do inexactness in the linear solver
177 if self.inexact_linear_ratio:
178 self.lin_tol = res * self.inexact_linear_ratio
180 if res < self.newton_tol:
181 break
183 # assemble dg
184 dg = Id - factor * (self.A + 1.0 / eps2 * sp.diags((1.0 - (nu + 1) * u**nu), offsets=0))
186 # newton update: u1 = u0 - g/dg
187 # u -= spsolve(dg, g)
188 u -= cg(
189 dg, g, x0=z, rtol=self.lin_tol, maxiter=self.lin_maxiter, atol=0, callback=self.work_counters['linear']
190 )[0]
191 # increase iteration count
192 n += 1
193 # print(n, res)
195 self.work_counters['newton']()
197 # if n == self.newton_maxiter:
198 # raise ProblemError('Newton did not converge after %i iterations, error is %s' % (n, res))
200 me = self.dtype_u(self.init)
201 me[:] = u.reshape(self.nvars)
203 self.newton_ncalls += 1
204 self.newton_itercount += n
206 return me
208 def eval_f(self, u, t):
209 """
210 Routine to evaluate the right-hand side of the problem.
212 Parameters
213 ----------
214 u : dtype_u
215 Current values of the numerical solution.
216 t : float
217 Current time of the numerical solution is computed (not used here).
219 Returns
220 -------
221 f : dtype_f
222 The right-hand side of the problem.
223 """
224 f = self.dtype_f(self.init)
225 v = u.flatten()
226 f[:] = (self.A.dot(v) + 1.0 / self.eps**2 * v * (1.0 - v**self.nu)).reshape(self.nvars)
228 self.work_counters['rhs']()
229 return f
231 def u_exact(self, t, u_init=None, t_init=None):
232 r"""
233 Routine to compute the exact solution at time :math:`t`.
235 Parameters
236 ----------
237 t : float
238 Time of the exact solution.
240 Returns
241 -------
242 me : dtype_u
243 The exact solution.
244 """
245 me = self.dtype_u(self.init, val=0.0)
246 if t > 0:
248 def eval_rhs(t, u):
249 return self.eval_f(u.reshape(self.init[0]), t).flatten()
251 me[:] = self.generate_scipy_reference_solution(eval_rhs, t, u_init, t_init)
253 else:
254 X, Y = np.meshgrid(self.xvalues, self.xvalues)
255 r2 = X**2 + Y**2
256 me[:] = np.tanh((self.radius - np.sqrt(r2)) / (np.sqrt(2) * self.eps))
258 return me
261# noinspection PyUnusedLocal
262class allencahn_semiimplicit(allencahn_fullyimplicit):
263 r"""
264 This class implements the two-dimensional Allen-Cahn equation with periodic boundary conditions :math:`u \in [-1, 1]^2`
266 .. math::
267 \frac{\partial u}{\partial t} = \Delta u + \frac{1}{\varepsilon^2} u (1 - u^\nu)
269 for constant parameter :math:`\nu`. Initial condition are circles of the form
271 .. math::
272 u({\bf x}, 0) = \tanh\left(\frac{r - \sqrt{x_i^2 + y_j^2}}{\sqrt{2}\varepsilon}\right)
274 for :math:`i, j=0,..,N-1`, where :math:`N` is the number of spatial grid points. For time-stepping, the problem is
275 treated in a *semi-implicit* way, i.e., the linear system containing the Laplacian is solved by the conjugate gradients
276 method, and the system containing the rest of the right-hand side is only evaluated at each time.
277 """
279 dtype_f = imex_mesh
281 def eval_f(self, u, t):
282 """
283 Routine to evaluate the right-hand side of the problem.
285 Parameters
286 ----------
287 u : dtype_u
288 Current values of the numerical solution.
289 t : float
290 Current time of the numerical solution is computed (not used here).
292 Returns
293 -------
294 f : dtype_f
295 The right-hand side of the problem.
296 """
297 f = self.dtype_f(self.init)
298 v = u.flatten()
299 f.impl[:] = self.A.dot(v).reshape(self.nvars)
300 f.expl[:] = (1.0 / self.eps**2 * v * (1.0 - v**self.nu)).reshape(self.nvars)
302 self.work_counters['rhs']()
303 return f
305 def solve_system(self, rhs, factor, u0, t):
306 r"""
307 Simple linear solver for :math:`(I-factor\cdot A)\vec{u}=\vec{rhs}`.
309 Parameters
310 ----------
311 rhs : dtype_f
312 Right-hand side for the linear system.
313 factor : float
314 Abbrev. for the local stepsize (or any other factor required).
315 u0 : dtype_u
316 Initial guess for the iterative solver.
317 t : float
318 Current time (e.g. for time-dependent BCs).
320 Returns
321 -------
322 me : dtype_u
323 The solution as mesh.
324 """
326 class context:
327 num_iter = 0
329 def callback(xk):
330 context.num_iter += 1
331 self.work_counters['linear']()
332 return context.num_iter
334 me = self.dtype_u(self.init)
336 Id = sp.eye(self.nvars[0] * self.nvars[1])
338 me[:] = cg(
339 Id - factor * self.A,
340 rhs.flatten(),
341 x0=u0.flatten(),
342 rtol=self.lin_tol,
343 maxiter=self.lin_maxiter,
344 atol=0,
345 callback=callback,
346 )[0].reshape(self.nvars)
348 self.lin_ncalls += 1
349 self.lin_itercount += context.num_iter
351 return me
353 def u_exact(self, t, u_init=None, t_init=None):
354 """
355 Routine to compute the exact solution at time t.
357 Parameters
358 ----------
359 t : float
360 Time of the exact solution.
362 Returns
363 -------
364 me : dtype_u
365 The exact solution.
366 """
367 me = self.dtype_u(self.init, val=0.0)
368 if t > 0:
370 def eval_rhs(t, u):
371 f = self.eval_f(u.reshape(self.init[0]), t)
372 return (f.impl + f.expl).flatten()
374 me[:] = self.generate_scipy_reference_solution(eval_rhs, t, u_init, t_init)
375 else:
376 me[:] = super().u_exact(t, u_init, t_init)
377 return me
380# noinspection PyUnusedLocal
381class allencahn_semiimplicit_v2(allencahn_fullyimplicit):
382 r"""
383 This class implements the two-dimensional Allen-Cahn (AC) equation with periodic boundary conditions :math:`u \in [-1, 1]^2`
385 .. math::
386 \frac{\partial u}{\partial t} = \Delta u + \frac{1}{\varepsilon^2} u (1 - u^\nu)
388 for constant parameter :math:`\nu`. Initial condition are circles of the form
390 .. math::
391 u({\bf x}, 0) = \tanh\left(\frac{r - \sqrt{x_i^2 + y_j^2}}{\sqrt{2}\varepsilon}\right)
393 for :math:`i, j=0,..,N-1`, where :math:`N` is the number of spatial grid points. For time-stepping, a special AC-splitting
394 is used to get a *semi-implicit* treatment of the problem: The term :math:`\Delta u - \frac{1}{\varepsilon^2} u^{\nu + 1}`
395 is handled implicitly and the nonlinear system including this part will be solved by Newton. :math:`\frac{1}{\varepsilon^2} u`
396 is only evaluated at each time.
397 """
399 dtype_f = imex_mesh
401 def eval_f(self, u, t):
402 """
403 Routine to evaluate the right-hand side of the problem.
405 Parameters
406 ----------
407 u : dtype_u
408 Current values of the numerical solution.
409 t : float
410 Current time of the numerical solution is computed.
412 Returns
413 -------
414 f : dtype_f
415 The right-hand side of the problem.
416 """
417 f = self.dtype_f(self.init)
418 v = u.flatten()
419 f.impl[:] = (self.A.dot(v) - 1.0 / self.eps**2 * v ** (self.nu + 1)).reshape(self.nvars)
420 f.expl[:] = (1.0 / self.eps**2 * v).reshape(self.nvars)
422 return f
424 def solve_system(self, rhs, factor, u0, t):
425 """
426 Simple Newton solver.
428 Parameters
429 ----------
430 rhs : dtype_f
431 Right-hand side for the nonlinear system.
432 factor : float
433 Abbrev. for the node-to-node stepsize (or any other factor required).
434 u0 : dtype_u
435 Initial guess for the iterative solver.
436 t : float
437 Current time (required here for the BC).
439 Returns
440 -------
441 me : dtype_u
442 The solution as mesh.
443 """
445 u = self.dtype_u(u0).flatten()
446 z = self.dtype_u(self.init, val=0.0).flatten()
447 nu = self.nu
448 eps2 = self.eps**2
450 Id = sp.eye(self.nvars[0] * self.nvars[1])
452 # start newton iteration
453 n = 0
454 res = 99
455 while n < self.newton_maxiter:
456 # form the function g with g(u) = 0
457 g = u - factor * (self.A.dot(u) - 1.0 / eps2 * u ** (nu + 1)) - rhs.flatten()
459 # if g is close to 0, then we are done
460 # res = np.linalg.norm(g, np.inf)
461 res = np.linalg.norm(g, np.inf)
463 if res < self.newton_tol:
464 break
466 # assemble dg
467 dg = Id - factor * (self.A - 1.0 / eps2 * sp.diags(((nu + 1) * u**nu), offsets=0))
469 # newton update: u1 = u0 - g/dg
470 # u -= spsolve(dg, g)
471 u -= cg(dg, g, x0=z, rtol=self.lin_tol, atol=0)[0]
472 # increase iteration count
473 n += 1
474 # print(n, res)
476 # if n == self.newton_maxiter:
477 # raise ProblemError('Newton did not converge after %i iterations, error is %s' % (n, res))
479 me = self.dtype_u(self.init)
480 me[:] = u.reshape(self.nvars)
482 self.newton_ncalls += 1
483 self.newton_itercount += n
485 return me
488# noinspection PyUnusedLocal
489class allencahn_multiimplicit(allencahn_fullyimplicit):
490 r"""
491 Example implementing the two-dimensional Allen-Cahn equation with periodic boundary conditions :math:`u \in [-1, 1]^2`
493 .. math::
494 \frac{\partial u}{\partial t} = \Delta u + \frac{1}{\varepsilon^2} u (1 - u^\nu)
496 for constant parameter :math:`\nu`. Initial condition are circles of the form
498 .. math::
499 u({\bf x}, 0) = \tanh\left(\frac{r - \sqrt{x_i^2 + y_j^2}}{\sqrt{2}\varepsilon}\right)
501 for :math:`i, j=0,..,N-1`, where :math:`N` is the number of spatial grid points. For time-stepping, the problem is
502 treated in *multi-implicit* fashion, i.e., the linear system containing the Laplacian is solved by the conjugate gradients
503 method, and the system containing the rest of the right-hand side will be solved by Newton's method.
504 """
506 dtype_f = comp2_mesh
508 def eval_f(self, u, t):
509 """
510 Routine to evaluate the right-hand side of the problem.
512 Parameters
513 ----------
514 u : dtype_u
515 Current values of the numerical solution.
516 t : float
517 Current time of the numerical solution is computed.
519 Returns
520 -------
521 f : dtype_f
522 The right-hand side of the problem.
523 """
524 f = self.dtype_f(self.init)
525 v = u.flatten()
526 f.comp1[:] = self.A.dot(v).reshape(self.nvars)
527 f.comp2[:] = (1.0 / self.eps**2 * v * (1.0 - v**self.nu)).reshape(self.nvars)
529 return f
531 def solve_system_1(self, rhs, factor, u0, t):
532 r"""
533 Simple linear solver for :math:`(I-factor\cdot A)\vec{u}=\vec{rhs}`.
535 Parameters
536 ----------
537 rhs : dtype_f
538 Right-hand side for the linear system.
539 factor : float
540 Abbrev. for the local stepsize (or any other factor required).
541 u0 : dtype_u
542 Initial guess for the iterative solver.
543 t : float
544 Current time (e.g. for time-dependent BCs).
546 Returns
547 -------
548 me : dtype_u
549 The solution as mesh.
550 """
552 class context:
553 num_iter = 0
555 def callback(xk):
556 context.num_iter += 1
557 return context.num_iter
559 me = self.dtype_u(self.init)
561 Id = sp.eye(self.nvars[0] * self.nvars[1])
563 me[:] = cg(
564 Id - factor * self.A,
565 rhs.flatten(),
566 x0=u0.flatten(),
567 rtol=self.lin_tol,
568 maxiter=self.lin_maxiter,
569 atol=0,
570 callback=callback,
571 )[0].reshape(self.nvars)
573 self.lin_ncalls += 1
574 self.lin_itercount += context.num_iter
576 return me
578 def solve_system_2(self, rhs, factor, u0, t):
579 """
580 Simple Newton solver.
582 Parameters
583 ----------
584 rhs : dtype_f
585 Right-hand side for the nonlinear system.
586 factor : float
587 Abbrev. for the node-to-node stepsize (or any other factor required).
588 u0 : dtype_u
589 Initial guess for the iterative solver.
590 t : float
591 Current time (required here for the BC).
593 Returns
594 -------
595 me : dtype_u
596 The solution as mesh.
597 """
599 u = self.dtype_u(u0).flatten()
600 z = self.dtype_u(self.init, val=0.0).flatten()
601 nu = self.nu
602 eps2 = self.eps**2
604 Id = sp.eye(self.nvars[0] * self.nvars[1])
606 # start newton iteration
607 n = 0
608 res = 99
609 while n < self.newton_maxiter:
610 # form the function g with g(u) = 0
611 g = u - factor * (1.0 / eps2 * u * (1.0 - u**nu)) - rhs.flatten()
613 # if g is close to 0, then we are done
614 res = np.linalg.norm(g, np.inf)
616 if res < self.newton_tol:
617 break
619 # assemble dg
620 dg = Id - factor * (1.0 / eps2 * sp.diags((1.0 - (nu + 1) * u**nu), offsets=0))
622 # newton update: u1 = u0 - g/dg
623 # u -= spsolve(dg, g)
624 u -= cg(dg, g, x0=z, rtol=self.lin_tol, atol=0)[0]
625 # increase iteration count
626 n += 1
627 # print(n, res)
629 # if n == self.newton_maxiter:
630 # raise ProblemError('Newton did not converge after %i iterations, error is %s' % (n, res))
632 me = self.dtype_u(self.init)
633 me[:] = u.reshape(self.nvars)
635 self.newton_ncalls += 1
636 self.newton_itercount += n
638 return me
641# noinspection PyUnusedLocal
642class allencahn_multiimplicit_v2(allencahn_fullyimplicit):
643 r"""
644 This class implements the two-dimensional Allen-Cahn (AC) equation with periodic boundary conditions :math:`u \in [-1, 1]^2`
646 .. math::
647 \frac{\partial u}{\partial t} = \Delta u + \frac{1}{\varepsilon^2} u (1 - u^\nu)
649 for constant parameter :math:`\nu`. The initial condition has the form of circles
651 .. math::
652 u({\bf x}, 0) = \tanh\left(\frac{r - \sqrt{x_i^2 + y_j^2}}{\sqrt{2}\varepsilon}\right)
654 for :math:`i, j=0,..,N-1`, where :math:`N` is the number of spatial grid points. For time-stepping, a special AC-splitting
655 is used here to get another kind of *semi-implicit* treatment of the problem: The term :math:`\Delta u - \frac{1}{\varepsilon^2} u^{\nu + 1}`
656 is handled implicitly and the nonlinear system including this part will be solved by Newton. :math:`\frac{1}{\varepsilon^2} u`
657 is solved by a linear solver provided by a ``SciPy`` routine.
658 """
660 dtype_f = comp2_mesh
662 def eval_f(self, u, t):
663 """
664 Routine to evaluate the right-hand side of the problem.
666 Parameters
667 ----------
668 u : dtype_u
669 Current values of the numerical solution.
670 t : float
671 Current time of the numerical solution is computed.
673 Returns
674 -------
675 f : dtype_f
676 The right-hand side of the problem.
677 """
678 f = self.dtype_f(self.init)
679 v = u.flatten()
680 f.comp1[:] = (self.A.dot(v) - 1.0 / self.eps**2 * v ** (self.nu + 1)).reshape(self.nvars)
681 f.comp2[:] = (1.0 / self.eps**2 * v).reshape(self.nvars)
683 return f
685 def solve_system_1(self, rhs, factor, u0, t):
686 """
687 Simple Newton solver.
689 Parameters
690 ----------
691 rhs : dtype_f
692 Right-hand side for the nonlinear system.
693 factor : float
694 Abbrev. for the node-to-node stepsize (or any other factor required).
695 u0 : dtype_u
696 Initial guess for the iterative solver.
697 t : float
698 Current time (required here for the BC).
700 Returns
701 ------
702 me : dtype_u
703 The solution as mesh.
704 """
706 u = self.dtype_u(u0).flatten()
707 z = self.dtype_u(self.init, val=0.0).flatten()
708 nu = self.nu
709 eps2 = self.eps**2
711 Id = sp.eye(self.nvars[0] * self.nvars[1])
713 # start newton iteration
714 n = 0
715 res = 99
716 while n < self.newton_maxiter:
717 # form the function g with g(u) = 0
718 g = u - factor * (self.A.dot(u) - 1.0 / eps2 * u ** (nu + 1)) - rhs.flatten()
720 # if g is close to 0, then we are done
721 res = np.linalg.norm(g, np.inf)
723 if res < self.newton_tol:
724 break
726 # assemble dg
727 dg = Id - factor * (self.A - 1.0 / eps2 * sp.diags(((nu + 1) * u**nu), offsets=0))
729 # newton update: u1 = u0 - g/dg
730 # u -= spsolve(dg, g)
731 u -= cg(
732 dg,
733 g,
734 x0=z,
735 rtol=self.lin_tol,
736 atol=0,
737 )[0]
738 # increase iteration count
739 n += 1
740 # print(n, res)
742 # if n == self.newton_maxiter:
743 # raise ProblemError('Newton did not converge after %i iterations, error is %s' % (n, res))
745 me = self.dtype_u(self.init)
746 me[:] = u.reshape(self.nvars)
748 self.newton_ncalls += 1
749 self.newton_itercount += n
751 return me
753 def solve_system_2(self, rhs, factor, u0, t):
754 r"""
755 Simple linear solver for :math:`(I-factor\cdot A)\vec{u}=\vec{rhs}`.
757 Parameters
758 ----------
759 rhs : dtype_f
760 Right-hand side for the linear system.
761 factor : float
762 Abbrev. for the local stepsize (or any other factor required).
763 u0 : dtype_u
764 Initial guess for the iterative solver.
765 t : float
766 Current time (e.g. for time-dependent BCs).
768 Returns
769 -------
770 me : dtype_u
771 The solution as mesh.
772 """
774 me = self.dtype_u(self.init)
776 me[:] = (1.0 / (1.0 - factor * 1.0 / self.eps**2) * rhs).reshape(self.nvars)
777 return me