Coverage for pySDC / implementations / problem_classes / AllenCahn_2D_FD.py: 96%
201 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-03-27 07:06 +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
10# http://www.personal.psu.edu/qud2/Res/Pre/dz09sisc.pdf
13# noinspection PyUnusedLocal
14class allencahn_fullyimplicit(Problem):
15 r"""
16 Example implementing the two-dimensional Allen-Cahn equation with periodic boundary conditions :math:`u \in [-1, 1]^2`
18 .. math::
19 \frac{\partial u}{\partial t} = \Delta u + \frac{1}{\varepsilon^2} u (1 - u^\nu)
21 for constant parameter :math:`\nu`. Initial condition are circles of the form
23 .. math::
24 u({\bf x}, 0) = \tanh\left(\frac{r - \sqrt{x_i^2 + y_j^2}}{\sqrt{2}\varepsilon}\right)
26 for :math:`i, j=0,..,N-1`, where :math:`N` is the number of spatial grid points. For time-stepping, the problem is
27 treated *fully-implicitly*, i.e., the nonlinear system is solved by Newton.
29 Parameters
30 ----------
31 nvars : tuple of int, optional
32 Number of unknowns in the problem, e.g. ``nvars=(128, 128)``.
33 nu : float, optional
34 Problem parameter :math:`\nu`.
35 eps : float, optional
36 Scaling parameter :math:`\varepsilon`.
37 newton_maxiter : int, optional
38 Maximum number of iterations for the Newton solver.
39 newton_tol : float, optional
40 Tolerance for Newton's method to terminate.
41 lin_tol : float, optional
42 Tolerance for linear solver to terminate.
43 lin_maxiter : int, optional
44 Maximum number of iterations for the linear solver.
45 radius : float, optional
46 Radius of the circles.
47 order : int, optional
48 Order of the finite difference matrix.
50 Attributes
51 ----------
52 A : scipy.spdiags
53 Second-order FD discretization of the 2D laplace operator.
54 dx : float
55 Distance between two spatial nodes (same for both directions).
56 xvalues : np.1darray
57 Spatial grid points, here both dimensions have the same grid points.
58 newton_itercount : int
59 Number of iterations of Newton solver.
60 lin_itercount
61 Number of iterations of linear solver.
62 newton_ncalls : int
63 Number of calls of Newton solver.
64 lin_ncalls : int
65 Number of calls of linear solver.
66 """
68 dtype_u = mesh
69 dtype_f = mesh
71 def __init__(
72 self,
73 nvars=(128, 128),
74 nu=2,
75 eps=0.04,
76 newton_maxiter=200,
77 newton_tol=1e-12,
78 lin_tol=1e-8,
79 lin_maxiter=100,
80 inexact_linear_ratio=None,
81 radius=0.25,
82 order=2,
83 ):
84 """Initialization routine"""
85 # we assert that nvars looks very particular here.. this will be necessary for coarsening in space later on
86 if len(nvars) != 2:
87 raise ProblemError('this is a 2d example, got %s' % nvars)
88 if nvars[0] != nvars[1]:
89 raise ProblemError('need a square domain, got %s' % nvars)
90 if nvars[0] % 2 != 0:
91 raise ProblemError('the setup requires nvars = 2^p per dimension')
93 # invoke super init, passing number of dofs, dtype_u and dtype_f
94 super().__init__((nvars, None, np.dtype('float64')))
95 self._makeAttributeAndRegister(
96 'nvars',
97 'nu',
98 'eps',
99 'radius',
100 'order',
101 localVars=locals(),
102 readOnly=True,
103 )
104 self._makeAttributeAndRegister(
105 'newton_maxiter',
106 'newton_tol',
107 'lin_tol',
108 'lin_maxiter',
109 'inexact_linear_ratio',
110 localVars=locals(),
111 readOnly=False,
112 )
114 # compute dx and get discretization matrix A
115 self.dx = 1.0 / self.nvars[0]
116 self.A, _ = problem_helper.get_finite_difference_matrix(
117 derivative=2,
118 order=self.order,
119 stencil_type='center',
120 dx=self.dx,
121 size=self.nvars[0],
122 dim=2,
123 bc='periodic',
124 )
125 self.xvalues = np.array([i * self.dx - 0.5 for i in range(self.nvars[0])])
127 self.newton_itercount = 0
128 self.lin_itercount = 0
129 self.newton_ncalls = 0
130 self.lin_ncalls = 0
132 self.work_counters['newton'] = WorkCounter()
133 self.work_counters['rhs'] = WorkCounter()
134 self.work_counters['linear'] = WorkCounter()
136 # noinspection PyTypeChecker
137 def solve_system(self, rhs, factor, u0, t):
138 """
139 Simple Newton solver.
141 Parameters
142 ----------
143 rhs : dtype_f
144 Right-hand side for the nonlinear system
145 factor : float
146 Abbrev. for the node-to-node stepsize (or any other factor required).
147 u0 : dtype_u
148 Initial guess for the iterative solver.
149 t : float
150 Current time (required here for the BC).
152 Returns
153 -------
154 me : dtype_u
155 The solution as mesh.
156 """
158 u = self.dtype_u(u0).flatten()
159 z = self.dtype_u(self.init, val=0.0).flatten()
160 nu = self.nu
161 eps2 = self.eps**2
163 Id = sp.eye(self.nvars[0] * self.nvars[1])
165 # start newton iteration
166 n = 0
167 res = 99
168 while n < self.newton_maxiter:
169 # form the function g with g(u) = 0
170 g = u - factor * (self.A.dot(u) + 1.0 / eps2 * u * (1.0 - u**nu)) - rhs.flatten()
172 # if g is close to 0, then we are done
173 res = np.linalg.norm(g, np.inf)
175 # do inexactness in the linear solver
176 if self.inexact_linear_ratio:
177 self.lin_tol = res * self.inexact_linear_ratio
179 if res < self.newton_tol:
180 break
182 # assemble dg
183 dg = Id - factor * (self.A + 1.0 / eps2 * sp.diags((1.0 - (nu + 1) * u**nu), offsets=0))
185 # newton update: u1 = u0 - g/dg
186 # u -= spsolve(dg, g)
187 u -= cg(
188 dg, g, x0=z, rtol=self.lin_tol, maxiter=self.lin_maxiter, atol=0, callback=self.work_counters['linear']
189 )[0]
190 # increase iteration count
191 n += 1
192 # print(n, res)
194 self.work_counters['newton']()
196 # if n == self.newton_maxiter:
197 # raise ProblemError('Newton did not converge after %i iterations, error is %s' % (n, res))
199 me = self.dtype_u(self.init)
200 me[:] = u.reshape(self.nvars)
202 self.newton_ncalls += 1
203 self.newton_itercount += n
205 return me
207 def eval_f(self, u, t):
208 """
209 Routine to evaluate the right-hand side of the problem.
211 Parameters
212 ----------
213 u : dtype_u
214 Current values of the numerical solution.
215 t : float
216 Current time of the numerical solution is computed (not used here).
218 Returns
219 -------
220 f : dtype_f
221 The right-hand side of the problem.
222 """
223 f = self.dtype_f(self.init)
224 v = u.flatten()
225 f[:] = (self.A.dot(v) + 1.0 / self.eps**2 * v * (1.0 - v**self.nu)).reshape(self.nvars)
227 self.work_counters['rhs']()
228 return f
230 def u_exact(self, t, u_init=None, t_init=None):
231 r"""
232 Routine to compute the exact solution at time :math:`t`.
234 Parameters
235 ----------
236 t : float
237 Time of the exact solution.
239 Returns
240 -------
241 me : dtype_u
242 The exact solution.
243 """
244 me = self.dtype_u(self.init, val=0.0)
245 if t > 0:
247 def eval_rhs(t, u):
248 return self.eval_f(u.reshape(self.init[0]), t).flatten()
250 me[:] = self.generate_scipy_reference_solution(eval_rhs, t, u_init, t_init)
252 else:
253 X, Y = np.meshgrid(self.xvalues, self.xvalues)
254 r2 = X**2 + Y**2
255 me[:] = np.tanh((self.radius - np.sqrt(r2)) / (np.sqrt(2) * self.eps))
257 return me
260# noinspection PyUnusedLocal
261class allencahn_semiimplicit(allencahn_fullyimplicit):
262 r"""
263 This class implements the two-dimensional Allen-Cahn equation with periodic boundary conditions :math:`u \in [-1, 1]^2`
265 .. math::
266 \frac{\partial u}{\partial t} = \Delta u + \frac{1}{\varepsilon^2} u (1 - u^\nu)
268 for constant parameter :math:`\nu`. Initial condition are circles of the form
270 .. math::
271 u({\bf x}, 0) = \tanh\left(\frac{r - \sqrt{x_i^2 + y_j^2}}{\sqrt{2}\varepsilon}\right)
273 for :math:`i, j=0,..,N-1`, where :math:`N` is the number of spatial grid points. For time-stepping, the problem is
274 treated in a *semi-implicit* way, i.e., the linear system containing the Laplacian is solved by the conjugate gradients
275 method, and the system containing the rest of the right-hand side is only evaluated at each time.
276 """
278 dtype_f = imex_mesh
280 def eval_f(self, u, t):
281 """
282 Routine to evaluate the right-hand side of the problem.
284 Parameters
285 ----------
286 u : dtype_u
287 Current values of the numerical solution.
288 t : float
289 Current time of the numerical solution is computed (not used here).
291 Returns
292 -------
293 f : dtype_f
294 The right-hand side of the problem.
295 """
296 f = self.dtype_f(self.init)
297 v = u.flatten()
298 f.impl[:] = self.A.dot(v).reshape(self.nvars)
299 f.expl[:] = (1.0 / self.eps**2 * v * (1.0 - v**self.nu)).reshape(self.nvars)
301 self.work_counters['rhs']()
302 return f
304 def solve_system(self, rhs, factor, u0, t):
305 r"""
306 Simple linear solver for :math:`(I-factor\cdot A)\vec{u}=\vec{rhs}`.
308 Parameters
309 ----------
310 rhs : dtype_f
311 Right-hand side for the linear system.
312 factor : float
313 Abbrev. for the local stepsize (or any other factor required).
314 u0 : dtype_u
315 Initial guess for the iterative solver.
316 t : float
317 Current time (e.g. for time-dependent BCs).
319 Returns
320 -------
321 me : dtype_u
322 The solution as mesh.
323 """
325 class context:
326 num_iter = 0
328 def callback(xk):
329 context.num_iter += 1
330 self.work_counters['linear']()
331 return context.num_iter
333 me = self.dtype_u(self.init)
335 Id = sp.eye(self.nvars[0] * self.nvars[1])
337 me[:] = cg(
338 Id - factor * self.A,
339 rhs.flatten(),
340 x0=u0.flatten(),
341 rtol=self.lin_tol,
342 maxiter=self.lin_maxiter,
343 atol=0,
344 callback=callback,
345 )[0].reshape(self.nvars)
347 self.lin_ncalls += 1
348 self.lin_itercount += context.num_iter
350 return me
352 def u_exact(self, t, u_init=None, t_init=None):
353 """
354 Routine to compute the exact solution at time t.
356 Parameters
357 ----------
358 t : float
359 Time of the exact solution.
361 Returns
362 -------
363 me : dtype_u
364 The exact solution.
365 """
366 me = self.dtype_u(self.init, val=0.0)
367 if t > 0:
369 def eval_rhs(t, u):
370 f = self.eval_f(u.reshape(self.init[0]), t)
371 return (f.impl + f.expl).flatten()
373 me[:] = self.generate_scipy_reference_solution(eval_rhs, t, u_init, t_init)
374 else:
375 me[:] = super().u_exact(t, u_init, t_init)
376 return me
379# noinspection PyUnusedLocal
380class allencahn_semiimplicit_v2(allencahn_fullyimplicit):
381 r"""
382 This class implements the two-dimensional Allen-Cahn (AC) equation with periodic boundary conditions :math:`u \in [-1, 1]^2`
384 .. math::
385 \frac{\partial u}{\partial t} = \Delta u + \frac{1}{\varepsilon^2} u (1 - u^\nu)
387 for constant parameter :math:`\nu`. Initial condition are circles of the form
389 .. math::
390 u({\bf x}, 0) = \tanh\left(\frac{r - \sqrt{x_i^2 + y_j^2}}{\sqrt{2}\varepsilon}\right)
392 for :math:`i, j=0,..,N-1`, where :math:`N` is the number of spatial grid points. For time-stepping, a special AC-splitting
393 is used to get a *semi-implicit* treatment of the problem: The term :math:`\Delta u - \frac{1}{\varepsilon^2} u^{\nu + 1}`
394 is handled implicitly and the nonlinear system including this part will be solved by Newton. :math:`\frac{1}{\varepsilon^2} u`
395 is only evaluated at each time.
396 """
398 dtype_f = imex_mesh
400 def eval_f(self, u, t):
401 """
402 Routine to evaluate the right-hand side of the problem.
404 Parameters
405 ----------
406 u : dtype_u
407 Current values of the numerical solution.
408 t : float
409 Current time of the numerical solution is computed.
411 Returns
412 -------
413 f : dtype_f
414 The right-hand side of the problem.
415 """
416 f = self.dtype_f(self.init)
417 v = u.flatten()
418 f.impl[:] = (self.A.dot(v) - 1.0 / self.eps**2 * v ** (self.nu + 1)).reshape(self.nvars)
419 f.expl[:] = (1.0 / self.eps**2 * v).reshape(self.nvars)
421 return f
423 def solve_system(self, rhs, factor, u0, t):
424 """
425 Simple Newton solver.
427 Parameters
428 ----------
429 rhs : dtype_f
430 Right-hand side for the nonlinear system.
431 factor : float
432 Abbrev. for the node-to-node stepsize (or any other factor required).
433 u0 : dtype_u
434 Initial guess for the iterative solver.
435 t : float
436 Current time (required here for the BC).
438 Returns
439 -------
440 me : dtype_u
441 The solution as mesh.
442 """
444 u = self.dtype_u(u0).flatten()
445 z = self.dtype_u(self.init, val=0.0).flatten()
446 nu = self.nu
447 eps2 = self.eps**2
449 Id = sp.eye(self.nvars[0] * self.nvars[1])
451 # start newton iteration
452 n = 0
453 res = 99
454 while n < self.newton_maxiter:
455 # form the function g with g(u) = 0
456 g = u - factor * (self.A.dot(u) - 1.0 / eps2 * u ** (nu + 1)) - rhs.flatten()
458 # if g is close to 0, then we are done
459 # res = np.linalg.norm(g, np.inf)
460 res = np.linalg.norm(g, np.inf)
462 if res < self.newton_tol:
463 break
465 # assemble dg
466 dg = Id - factor * (self.A - 1.0 / eps2 * sp.diags(((nu + 1) * u**nu), offsets=0))
468 # newton update: u1 = u0 - g/dg
469 # u -= spsolve(dg, g)
470 u -= cg(dg, g, x0=z, rtol=self.lin_tol, atol=0)[0]
471 # increase iteration count
472 n += 1
473 # print(n, res)
475 # if n == self.newton_maxiter:
476 # raise ProblemError('Newton did not converge after %i iterations, error is %s' % (n, res))
478 me = self.dtype_u(self.init)
479 me[:] = u.reshape(self.nvars)
481 self.newton_ncalls += 1
482 self.newton_itercount += n
484 return me
487# noinspection PyUnusedLocal
488class allencahn_multiimplicit(allencahn_fullyimplicit):
489 r"""
490 Example implementing the two-dimensional Allen-Cahn equation with periodic boundary conditions :math:`u \in [-1, 1]^2`
492 .. math::
493 \frac{\partial u}{\partial t} = \Delta u + \frac{1}{\varepsilon^2} u (1 - u^\nu)
495 for constant parameter :math:`\nu`. Initial condition are circles of the form
497 .. math::
498 u({\bf x}, 0) = \tanh\left(\frac{r - \sqrt{x_i^2 + y_j^2}}{\sqrt{2}\varepsilon}\right)
500 for :math:`i, j=0,..,N-1`, where :math:`N` is the number of spatial grid points. For time-stepping, the problem is
501 treated in *multi-implicit* fashion, i.e., the linear system containing the Laplacian is solved by the conjugate gradients
502 method, and the system containing the rest of the right-hand side will be solved by Newton's method.
503 """
505 dtype_f = comp2_mesh
507 def eval_f(self, u, t):
508 """
509 Routine to evaluate the right-hand side of the problem.
511 Parameters
512 ----------
513 u : dtype_u
514 Current values of the numerical solution.
515 t : float
516 Current time of the numerical solution is computed.
518 Returns
519 -------
520 f : dtype_f
521 The right-hand side of the problem.
522 """
523 f = self.dtype_f(self.init)
524 v = u.flatten()
525 f.comp1[:] = self.A.dot(v).reshape(self.nvars)
526 f.comp2[:] = (1.0 / self.eps**2 * v * (1.0 - v**self.nu)).reshape(self.nvars)
528 return f
530 def solve_system_1(self, rhs, factor, u0, t):
531 r"""
532 Simple linear solver for :math:`(I-factor\cdot A)\vec{u}=\vec{rhs}`.
534 Parameters
535 ----------
536 rhs : dtype_f
537 Right-hand side for the linear system.
538 factor : float
539 Abbrev. for the local stepsize (or any other factor required).
540 u0 : dtype_u
541 Initial guess for the iterative solver.
542 t : float
543 Current time (e.g. for time-dependent BCs).
545 Returns
546 -------
547 me : dtype_u
548 The solution as mesh.
549 """
551 class context:
552 num_iter = 0
554 def callback(xk):
555 context.num_iter += 1
556 return context.num_iter
558 me = self.dtype_u(self.init)
560 Id = sp.eye(self.nvars[0] * self.nvars[1])
562 me[:] = cg(
563 Id - factor * self.A,
564 rhs.flatten(),
565 x0=u0.flatten(),
566 rtol=self.lin_tol,
567 maxiter=self.lin_maxiter,
568 atol=0,
569 callback=callback,
570 )[0].reshape(self.nvars)
572 self.lin_ncalls += 1
573 self.lin_itercount += context.num_iter
575 return me
577 def solve_system_2(self, rhs, factor, u0, t):
578 """
579 Simple Newton solver.
581 Parameters
582 ----------
583 rhs : dtype_f
584 Right-hand side for the nonlinear system.
585 factor : float
586 Abbrev. for the node-to-node stepsize (or any other factor required).
587 u0 : dtype_u
588 Initial guess for the iterative solver.
589 t : float
590 Current time (required here for the BC).
592 Returns
593 -------
594 me : dtype_u
595 The solution as mesh.
596 """
598 u = self.dtype_u(u0).flatten()
599 z = self.dtype_u(self.init, val=0.0).flatten()
600 nu = self.nu
601 eps2 = self.eps**2
603 Id = sp.eye(self.nvars[0] * self.nvars[1])
605 # start newton iteration
606 n = 0
607 res = 99
608 while n < self.newton_maxiter:
609 # form the function g with g(u) = 0
610 g = u - factor * (1.0 / eps2 * u * (1.0 - u**nu)) - rhs.flatten()
612 # if g is close to 0, then we are done
613 res = np.linalg.norm(g, np.inf)
615 if res < self.newton_tol:
616 break
618 # assemble dg
619 dg = Id - factor * (1.0 / eps2 * sp.diags((1.0 - (nu + 1) * u**nu), offsets=0))
621 # newton update: u1 = u0 - g/dg
622 # u -= spsolve(dg, g)
623 u -= cg(dg, g, x0=z, rtol=self.lin_tol, atol=0)[0]
624 # increase iteration count
625 n += 1
626 # print(n, res)
628 # if n == self.newton_maxiter:
629 # raise ProblemError('Newton did not converge after %i iterations, error is %s' % (n, res))
631 me = self.dtype_u(self.init)
632 me[:] = u.reshape(self.nvars)
634 self.newton_ncalls += 1
635 self.newton_itercount += n
637 return me
640# noinspection PyUnusedLocal
641class allencahn_multiimplicit_v2(allencahn_fullyimplicit):
642 r"""
643 This class implements the two-dimensional Allen-Cahn (AC) equation with periodic boundary conditions :math:`u \in [-1, 1]^2`
645 .. math::
646 \frac{\partial u}{\partial t} = \Delta u + \frac{1}{\varepsilon^2} u (1 - u^\nu)
648 for constant parameter :math:`\nu`. The initial condition has the form of circles
650 .. math::
651 u({\bf x}, 0) = \tanh\left(\frac{r - \sqrt{x_i^2 + y_j^2}}{\sqrt{2}\varepsilon}\right)
653 for :math:`i, j=0,..,N-1`, where :math:`N` is the number of spatial grid points. For time-stepping, a special AC-splitting
654 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}`
655 is handled implicitly and the nonlinear system including this part will be solved by Newton. :math:`\frac{1}{\varepsilon^2} u`
656 is solved by a linear solver provided by a ``SciPy`` routine.
657 """
659 dtype_f = comp2_mesh
661 def eval_f(self, u, t):
662 """
663 Routine to evaluate the right-hand side of the problem.
665 Parameters
666 ----------
667 u : dtype_u
668 Current values of the numerical solution.
669 t : float
670 Current time of the numerical solution is computed.
672 Returns
673 -------
674 f : dtype_f
675 The right-hand side of the problem.
676 """
677 f = self.dtype_f(self.init)
678 v = u.flatten()
679 f.comp1[:] = (self.A.dot(v) - 1.0 / self.eps**2 * v ** (self.nu + 1)).reshape(self.nvars)
680 f.comp2[:] = (1.0 / self.eps**2 * v).reshape(self.nvars)
682 return f
684 def solve_system_1(self, rhs, factor, u0, t):
685 """
686 Simple Newton solver.
688 Parameters
689 ----------
690 rhs : dtype_f
691 Right-hand side for the nonlinear system.
692 factor : float
693 Abbrev. for the node-to-node stepsize (or any other factor required).
694 u0 : dtype_u
695 Initial guess for the iterative solver.
696 t : float
697 Current time (required here for the BC).
699 Returns
700 ------
701 me : dtype_u
702 The solution as mesh.
703 """
705 u = self.dtype_u(u0).flatten()
706 z = self.dtype_u(self.init, val=0.0).flatten()
707 nu = self.nu
708 eps2 = self.eps**2
710 Id = sp.eye(self.nvars[0] * self.nvars[1])
712 # start newton iteration
713 n = 0
714 res = 99
715 while n < self.newton_maxiter:
716 # form the function g with g(u) = 0
717 g = u - factor * (self.A.dot(u) - 1.0 / eps2 * u ** (nu + 1)) - rhs.flatten()
719 # if g is close to 0, then we are done
720 res = np.linalg.norm(g, np.inf)
722 if res < self.newton_tol:
723 break
725 # assemble dg
726 dg = Id - factor * (self.A - 1.0 / eps2 * sp.diags(((nu + 1) * u**nu), offsets=0))
728 # newton update: u1 = u0 - g/dg
729 # u -= spsolve(dg, g)
730 u -= cg(
731 dg,
732 g,
733 x0=z,
734 rtol=self.lin_tol,
735 atol=0,
736 )[0]
737 # increase iteration count
738 n += 1
739 # print(n, res)
741 # if n == self.newton_maxiter:
742 # raise ProblemError('Newton did not converge after %i iterations, error is %s' % (n, res))
744 me = self.dtype_u(self.init)
745 me[:] = u.reshape(self.nvars)
747 self.newton_ncalls += 1
748 self.newton_itercount += n
750 return me
752 def solve_system_2(self, rhs, factor, u0, t):
753 r"""
754 Simple linear solver for :math:`(I-factor\cdot A)\vec{u}=\vec{rhs}`.
756 Parameters
757 ----------
758 rhs : dtype_f
759 Right-hand side for the linear system.
760 factor : float
761 Abbrev. for the local stepsize (or any other factor required).
762 u0 : dtype_u
763 Initial guess for the iterative solver.
764 t : float
765 Current time (e.g. for time-dependent BCs).
767 Returns
768 -------
769 me : dtype_u
770 The solution as mesh.
771 """
773 me = self.dtype_u(self.init)
775 me[:] = (1.0 / (1.0 - factor * 1.0 / self.eps**2) * rhs).reshape(self.nvars)
776 return me