import numpy as np
from petsc4py import PETSc
from pySDC.core.problem import Problem
from pySDC.implementations.datatype_classes.petsc_vec import petsc_vec, petsc_vec_imex, petsc_vec_comp2
[docs]
class Fisher_full(object):
r"""
Helper class to generate residual and Jacobian matrix for PETSc's nonlinear solver SNES.
Parameters
----------
da : DMDA object
Object of ``PETSc``.
prob : problem instance
Contains problem information for ``PETSc``.
factor : float
Temporal factor :math:`\Delta t Q_\Delta`.
dx : float
Grid spacing in x direction.
Attributes
----------
localX : PETSc vector object
Local vector for ``PETSc``.
xs, xe : int
Defines the ranges for spatial domain.
mx : tuple
Get sizes for the vector containing the spatial points.
row : PETSc matrix stencil object
Row for a matrix.
col : PETSc matrix stencil object
Column for a matrix.
"""
def __init__(self, da, prob, factor, dx):
"""Initialization routine"""
assert da.getDim() == 1
self.da = da
self.factor = factor
self.dx = dx
self.prob = prob
self.localX = da.createLocalVec()
self.xs, self.xe = self.da.getRanges()[0]
self.mx = self.da.getSizes()[0]
self.row = PETSc.Mat.Stencil()
self.col = PETSc.Mat.Stencil()
[docs]
class Fisher_reaction(object):
r"""
Helper class to generate residual and Jacobian matrix for ``PETSc``'s nonlinear solver SNES.
Parameters
----------
da : DMDA object
Object of ``PETSc``.
prob : problem instance
Contains problem information for ``PETSc``.
factor : float
Temporal factor :math:`\Delta t Q_\Delta`.
dx : float
Grid spacing in x direction.
Attributes
----------
localX : PETSc vector object
Local vector for ``PETSc``.
"""
def __init__(self, da, prob, factor):
"""Initialization routine"""
assert da.getDim() == 1
self.da = da
self.prob = prob
self.factor = factor
self.localX = da.createLocalVec()
[docs]
class petsc_fisher_multiimplicit(Problem):
r"""
The following one-dimensional problem is an example of a reaction-diffusion equation with traveling waves, and can
be seen as a generalized Fisher equation. This class implements a special case of the Kolmogorov-Petrovskii-Piskunov
problem [1]_ using periodic boundary conditions
.. math::
\frac{\partial u}{\partial t} = \Delta u + \lambda_0^2 u (1 - u^\nu)
with exact solution
.. math::
u(x, 0) = \left[
1 + \left(2^{\nu / 2} - 1\right) \exp\left(-(\nu / 2)\sigma_1 x + 2 \lambda_1 t\right)
\right]^{-2 / \nu}
for :math:`x \in \mathbb{R}`, and
.. math::
\sigma_1 = \lambda_1 - \sqrt{\lambda_1^2 - \lambda_0^2},
.. math::
\lambda_1 = \frac{\lambda_0}{2} \left[
\left(1 + \frac{\nu}{2}\right)^{1/2} + \left(1 + \frac{\nu}{2}\right)^{-1/2}
\right].
This class is implemented to be solved in spatial using ``PETSc`` [2]_, [3]_. For time-stepping, the problem will be solved
*multi-implicitly*.
Parameters
----------
nvars : int, optional
Spatial resolution.
lambda0 : float, optional
Problem parameter : math:`\lambda_0`.
nu : float, optional
Problem parameter :math:`\nu`.
interval : tuple, optional
Defines the spatial domain.
comm : PETSc.COMM_WORLD, optional
Communicator for PETSc.
lsol_tol : float, optional
Tolerance for linear solver to terminate.
nlsol_tol : float, optional
Tolerance for nonlinear solver to terminate.
lsol_maxiter : int, optional
Maximum number of iterations for linear solver.
nlsol_maxiter : int, optional
Maximum number of iterations for nonlinear solver.
Attributes
----------
dx : float
Mesh grid width.
xs, xe : int
Define the ranges.
A : PETSc matrix object
Discretization matrix.
localX : PETSc vector object
Local vector for solution.
ksp : PETSc solver object
Linear solver object.
snes : PETSc solver object
Nonlinear solver object.
F : PETSc vector object
Global vector.
J : PETSc matrix object
Jacobi matrix.
References
----------
.. [1] Z. Feng. Traveling wave behavior for a generalized fisher equation. Chaos Solitons Fractals 38(2),
481–488 (2008).
.. [2] PETSc Web page. Satish Balay et al. https://petsc.org/ (2023).
.. [3] Parallel distributed computing using Python. Lisandro D. Dalcin, Rodrigo R. Paz, Pablo A. Kler,
Alejandro Cosimo. Advances in Water Resources (2011).
"""
dtype_u = petsc_vec
dtype_f = petsc_vec_comp2
def __init__(
self,
nvars,
lambda0,
nu,
interval,
comm=PETSc.COMM_WORLD,
lsol_tol=1e-10,
nlsol_tol=1e-10,
lsol_maxiter=None,
nlsol_maxiter=None,
):
"""Initialization routine"""
# create DMDA object which will be used for all grid operations
da = PETSc.DMDA().create([nvars], dof=1, stencil_width=1, comm=comm)
# invoke super init, passing number of dofs, dtype_u and dtype_f
super().__init__(init=da)
self._makeAttributeAndRegister(
'nvars',
'lambda0',
'nu',
'interval',
'comm',
'lsol_tol',
'nlsol_tol',
'lsol_maxiter',
'nlsol_maxiter',
localVars=locals(),
readOnly=True,
)
# compute dx and get local ranges
self.dx = (self.interval[1] - self.interval[0]) / (self.nvars - 1)
(self.xs, self.xe) = self.init.getRanges()[0]
# compute discretization matrix A and identity
self.A = self.__get_A()
self.localX = self.init.createLocalVec()
# setup linear solver
self.ksp = PETSc.KSP()
self.ksp.create(comm=self.comm)
self.ksp.setType('cg')
pc = self.ksp.getPC()
pc.setType('ilu')
self.ksp.setInitialGuessNonzero(True)
self.ksp.setFromOptions()
self.ksp.setTolerances(rtol=self.lsol_tol, atol=self.lsol_tol, max_it=self.lsol_maxiter)
self.ksp_itercount = 0
self.ksp_ncalls = 0
# setup nonlinear solver
self.snes = PETSc.SNES()
self.snes.create(comm=self.comm)
if self.nlsol_maxiter <= 1:
self.snes.setType('ksponly')
self.snes.getKSP().setType('cg')
pc = self.snes.getKSP().getPC()
pc.setType('ilu')
# self.snes.setType('ngmres')
self.snes.setFromOptions()
self.snes.setTolerances(
rtol=self.nlsol_tol,
atol=self.nlsol_tol,
stol=self.nlsol_tol,
max_it=self.nlsol_maxiter,
)
self.snes_itercount = 0
self.snes_ncalls = 0
self.F = self.init.createGlobalVec()
self.J = self.init.createMatrix()
def __get_A(self):
r"""
Helper function to assemble ``PETSc`` matrix A.
Returns
-------
A : PETSc matrix object
Discretization matrix.
"""
# create matrix and set basic options
A = self.init.createMatrix()
A.setType('aij') # sparse
A.setFromOptions()
A.setPreallocationNNZ((3, 3))
A.setUp()
# fill matrix
A.zeroEntries()
row = PETSc.Mat.Stencil()
col = PETSc.Mat.Stencil()
mx = self.init.getSizes()[0]
(xs, xe) = self.init.getRanges()[0]
for i in range(xs, xe):
row.i = i
row.field = 0
if i == 0 or i == mx - 1:
A.setValueStencil(row, row, 1.0)
else:
diag = -2.0 / self.dx**2
for index, value in [
(i - 1, 1.0 / self.dx**2),
(i, diag),
(i + 1, 1.0 / self.dx**2),
]:
col.i = index
col.field = 0
A.setValueStencil(row, col, value)
A.assemble()
return A
[docs]
def get_sys_mat(self, factor):
"""
Helper function to assemble the system matrix of the linear problem.
Parameters
----------
factor : float
Factor to define the system matrix.
Returns
-------
A : PETSc matrix object
Matrix for the system to solve.
"""
# create matrix and set basic options
A = self.init.createMatrix()
A.setType('aij') # sparse
A.setFromOptions()
A.setPreallocationNNZ((3, 3))
A.setUp()
# fill matrix
A.zeroEntries()
row = PETSc.Mat.Stencil()
col = PETSc.Mat.Stencil()
mx = self.init.getSizes()[0]
(xs, xe) = self.init.getRanges()[0]
for i in range(xs, xe):
row.i = i
row.field = 0
if i == 0 or i == mx - 1:
A.setValueStencil(row, row, 1.0)
else:
diag = 1.0 + factor * 2.0 / self.dx**2
for index, value in [
(i - 1, -factor / self.dx**2),
(i, diag),
(i + 1, -factor / self.dx**2),
]:
col.i = index
col.field = 0
A.setValueStencil(row, col, value)
A.assemble()
return A
[docs]
def eval_f(self, u, t):
"""
Routine to evaluate the right-hand side of the problem.
Parameters
----------
u : dtype_u
Current values of the numerical solution.
t : float
Current time the numerical solution is computed.
Returns
-------
f : dtype_f
Right-hand side of the problem.
"""
f = self.dtype_f(self.init)
self.A.mult(u, f.comp1)
fa1 = self.init.getVecArray(f.comp1)
fa1[0] = 0
fa1[-1] = 0
fa2 = self.init.getVecArray(f.comp2)
xa = self.init.getVecArray(u)
for i in range(self.xs, self.xe):
fa2[i] = self.lambda0**2 * xa[i] * (1 - xa[i] ** self.nu)
fa2[0] = 0
fa2[-1] = 0
return f
[docs]
def solve_system_1(self, rhs, factor, u0, t):
r"""
Linear solver for :math:`(I - factor \cdot A)\vec{u} = \vec{rhs}`.
Parameters
----------
rhs : dtype_f
Right-hand side for the linear system.
factor : float
Abbrev. for the local stepsize (or any other factor required).
u0 : dtype_u
Initial guess for the iterative solver.
t : float
Current time (e.g. for time-dependent BCs).
Returns
-------
me : dtype_u
Solution as mesh.
"""
me = self.dtype_u(u0)
self.ksp.setOperators(self.get_sys_mat(factor))
self.ksp.solve(rhs, me)
self.ksp_itercount += self.ksp.getIterationNumber()
self.ksp_ncalls += 1
return me
[docs]
def solve_system_2(self, rhs, factor, u0, t):
r"""
Nonlinear solver for :math:`(I - factor \cdot F)(\vec{u}) = \vec{rhs}`.
Parameters
----------
rhs : dtype_f
Right-hand side for the linear system.
factor : float
Abbrev. for the local stepsize (or any other factor required).
u0 : dtype_u
Initial guess for the iterative solver.
t : float
Current time (e.g. for time-dependent BCs).
Returns
-------
me : dtype_u
Solution as mesh.
"""
me = self.dtype_u(u0)
target = Fisher_reaction(self.init, self, factor)
# assign residual function and Jacobian
F = self.init.createGlobalVec()
self.snes.setFunction(target.formFunction, F)
J = self.init.createMatrix()
self.snes.setJacobian(target.formJacobian, J)
self.snes.solve(rhs, me)
self.snes_itercount += self.snes.getIterationNumber()
self.snes_ncalls += 1
return me
[docs]
def u_exact(self, t):
r"""
Routine to compute the exact solution at time :math:`t`.
Parameters
----------
t : float
Time of the exact solution.
Returns
-------
me : dtype_u
Exact solution.
"""
lam1 = self.lambda0 / 2.0 * ((self.nu / 2.0 + 1) ** 0.5 + (self.nu / 2.0 + 1) ** (-0.5))
sig1 = lam1 - np.sqrt(lam1**2 - self.lambda0**2)
me = self.dtype_u(self.init)
xa = self.init.getVecArray(me)
for i in range(self.xs, self.xe):
xa[i] = (
1
+ (2 ** (self.nu / 2.0) - 1)
* np.exp(-self.nu / 2.0 * sig1 * (self.interval[0] + (i + 1) * self.dx + 2 * lam1 * t))
) ** (-2.0 / self.nu)
return me
[docs]
class petsc_fisher_fullyimplicit(petsc_fisher_multiimplicit):
r"""
The following one-dimensional problem is an example of a reaction-diffusion equation with traveling waves, and can
be seen as a generalized Fisher equation. This class implements a special case of the Kolmogorov-Petrovskii-Piskunov
problem [1]_ using periodic boundary conditions
.. math::
\frac{\partial u}{\partial t} = \Delta u + \lambda_0^2 u (1 - u^\nu)
with exact solution
.. math::
u(x, 0) = \left[
1 + \left(2^{\nu / 2} - 1\right) \exp\left(-(\nu / 2)\sigma_1 x + 2 \lambda_1 t\right)
\right]^{-2 / \nu}
for :math:`x \in \mathbb{R}`, and
.. math::
\sigma_1 = \lambda_1 - \sqrt{\lambda_1^2 - \lambda_0^2},
.. math::
\lambda_1 = \frac{\lambda_0}{2} \left[
\left(1 + \frac{\nu}{2}\right)^{1/2} + \left(1 + \frac{\nu}{2}\right)^{-1/2}
\right].
This class is implemented to be solved in spatial using ``PETSc`` [2]_, [3]_. For time-stepping, the problem is treated
*fully-implicitly*.
"""
dtype_f = petsc_vec
[docs]
def eval_f(self, u, t):
"""
Routine to evaluate the right-hand side of the problem.
Parameters
----------
u : dtype_u
Current values of the numerical solution.
t : float
Current time the numerical solution is computed.
Returns
-------
f : dtype_f
Right-hand side of the problem.
"""
f = self.dtype_f(self.init)
self.A.mult(u, f)
fa2 = self.init.getVecArray(f)
xa = self.init.getVecArray(u)
for i in range(self.xs, self.xe):
fa2[i] += self.lambda0**2 * xa[i] * (1 - xa[i] ** self.nu)
fa2[0] = 0
fa2[-1] = 0
return f
[docs]
def solve_system(self, rhs, factor, u0, t):
r"""
Nonlinear solver for :math:`(I - factor \cdot F)(\vec{u}) = \vec{rhs}`.
Parameters
----------
rhs : dtype_f
Right-hand side for the linear system.
factor : float
Abbrev. for the local stepsize (or any other factor required).
u0 : dtype_u
Initial guess for the iterative solver.
t : float
Current time (e.g. for time-dependent BCs).
Returns
-------
me : dtype_u
Solution as mesh.
"""
me = self.dtype_u(u0)
target = Fisher_full(self.init, self, factor, self.dx)
# assign residual function and Jacobian
self.snes.setFunction(target.formFunction, self.F)
self.snes.setJacobian(target.formJacobian, self.J)
self.snes.solve(rhs, me)
self.snes_itercount += self.snes.getIterationNumber()
self.snes_ncalls += 1
return me
[docs]
class petsc_fisher_semiimplicit(petsc_fisher_multiimplicit):
r"""
The following one-dimensional problem is an example of a reaction-diffusion equation with traveling waves, and can
be seen as a generalized Fisher equation. This class implements a special case of the Kolmogorov-Petrovskii-Piskunov
problem [1]_ using periodic boundary conditions
.. math::
\frac{\partial u}{\partial t} = \Delta u + \lambda_0^2 u (1 - u^\nu)
with exact solution
.. math::
u(x, 0) = \left[
1 + \left(2^{\nu / 2} - 1\right) \exp\left(-(\nu / 2)\sigma_1 x + 2 \lambda_1 t\right)
\right]^{-2 / \nu}
for :math:`x \in \mathbb{R}`, and
.. math::
\sigma_1 = \lambda_1 - \sqrt{\lambda_1^2 - \lambda_0^2},
.. math::
\lambda_1 = \frac{\lambda_0}{2} \left[
\left(1 + \frac{\nu}{2}\right)^{1/2} + \left(1 + \frac{\nu}{2}\right)^{-1/2}
\right].
This class is implemented to be solved in spatial using ``PETSc`` [2]_, [3]_. For time-stepping, the problem here will be
solved in a *semi-implicit* way.
"""
dtype_f = petsc_vec_imex
[docs]
def eval_f(self, u, t):
"""
Routine to evaluate the right-hand side of the problem.
Parameters
----------
u : dtype_u
Current values of the numerical solution.
t : float
Current time the numerical solution is computed.
Returns
-------
f : dtype_f
Right-hand side of the problem.
"""
f = self.dtype_f(self.init)
self.A.mult(u, f.impl)
fa1 = self.init.getVecArray(f.impl)
fa1[0] = 0
fa1[-1] = 0
fa2 = self.init.getVecArray(f.expl)
xa = self.init.getVecArray(u)
for i in range(self.xs, self.xe):
fa2[i] = self.lambda0**2 * xa[i] * (1 - xa[i] ** self.nu)
fa2[0] = 0
fa2[-1] = 0
return f
[docs]
def solve_system(self, rhs, factor, u0, t):
r"""
Simple linear solver for :math:`(I-factor \cdot A)\vec{u} = \vec{rhs}`.
Parameters
----------
rhs : dtype_f
Right-hand side for the linear system.
factor : float
Abbrev. for the local stepsize (or any other factor required).
u0 : dtype_u
Initial guess for the iterative solver.
t : float
Current time (e.g. for time-dependent BCs).
Returns
-------
me : dtype_u
Solution as mesh.
"""
me = self.dtype_u(u0)
self.ksp.setOperators(self.get_sys_mat(factor))
self.ksp.solve(rhs, me)
self.ksp_itercount += self.ksp.getIterationNumber()
self.ksp_ncalls += 1
return me