Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Added preconditioner based on Taylor expansion #158

Open
wants to merge 10 commits into
base: master
Choose a base branch
from
61 changes: 60 additions & 1 deletion pySDC/core/Sweeper.py
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@
from pySDC.implementations.collocation_classes.equidistant_right import EquidistantNoLeft
from pySDC.implementations.collocation_classes.gauss_lobatto import CollGaussLobatto
from pySDC.implementations.collocation_classes.gauss_radau_right import CollGaussRadau_Right
from pySDC.helpers.preconditioner_helper import get_linear_multistep_method


# short helper class to add params as attributes
Expand Down Expand Up @@ -200,8 +201,66 @@ def rho(x):
'implemented')
QDmat[1:, 1:] = np.diag(x)
self.parallelizable = True
elif qd_type == 'LMMEuler':
'''
Do a Taylor expansion using u0 and f(u_m) backwards in time from tau_m and solve for u_m.
We take as many values as we have available except for the rhs. of the initial conditions, meaning the order
of the Taylor expansion is m.
That means for the first node, we get backward Euler, but for later nodes we get higher order
solutions and in some cases we can get more than one order per sweep.
'''
for i in range(1, len(self.coll.nodes) + 1):
t_expand = self.coll.nodes[i - 1]
h = np.append([0], self.coll.nodes[:i]) - t_expand # time difference to where we expand about

u_signature = np.zeros_like(h)
u_signature[0] = 1

f_signature = np.ones_like(h)
f_signature[0] = 0

u_coeff, f_coeff = get_linear_multistep_method(h, u_signature, f_signature)

QDmat[i, 0: i + 1] = f_coeff
elif qd_type == 'LMM':
'''
Do a Taylor expansion using u0 and f(u_m) backwards in time from tau_m and solve for u_m.
We take as many values as we have available, meaning the order of the Taylor expansion is m + 1.
That means for the first node, we get trapezoidal rule, but for later nodes we get higher order
solutions and in some cases we can get more than one order per sweep.
'''
for i in range(1, len(self.coll.nodes) + 1):
t_expand = self.coll.nodes[i - 1]
h = np.append([0], self.coll.nodes[:i]) - t_expand # time difference to where we expand about

u_signature = np.zeros_like(h)
u_signature[0] = 1

f_signature = np.ones_like(h)

u_coeff, f_coeff = get_linear_multistep_method(h, u_signature, f_signature)

QDmat[i, 0: i + 1] = f_coeff
elif qd_type == 'LMMpar':
'''
Trapezoidal rule between initial conditions and the node you want to compute the solution at.
'''
for i in range(1, len(self.coll.nodes) + 1):
t_expand = self.coll.nodes[i - 1]
h = np.append([0], self.coll.nodes[:i]) - t_expand # time difference to where we expand about

u_signature = np.zeros_like(h)
u_signature[0] = 1

f_signature = np.zeros_like(h)
f_signature[-1] = 1
f_signature[0] = 1

u_coeff, f_coeff = get_linear_multistep_method(h, u_signature, f_signature)

QDmat[i, 0: i + 1] = f_coeff
else:
raise NotImplementedError('qd_type implicit not implemented')
raise NotImplementedError(f'qd_type implicit "{qd_type}" not implemented')
# check if we got not more than a lower triangular matrix
np.testing.assert_array_equal(np.triu(QDmat, k=1), np.zeros(QDmat.shape),
err_msg='Lower triangular matrix expected!')
Expand Down
220 changes: 220 additions & 0 deletions pySDC/helpers/preconditioner_helper.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,220 @@
import numpy as np
from scipy.special import factorial


def get_linear_multistep_method(steps, u_signature, f_signature, checks=False):
'''
Derive a general linear multistep method from step sizes and a signature. This function will provide a consistent
linear multistep method by cancelling terms in Taylor expansions, but note that this must not be stable!

The resulting coefficients must be multiplied with the corresponding value of u or f and then all must be summed to
get a numerical solution to the initial value problem.

Since we can cancel as many terms in the Taylor expansion as we have entries in the signature, that will also be
the order of consistency of our method.

We check if our method is consistent and zero stable, which together means it is convergent.
However, some of the methods that we generate are not A stable. As it turns out, according to Dahlquist's second
barrier theorem, there are no A-stable LMMs of order greater than 2.

Args:
steps (list): The step sizes between the multiple steps that are used
u_signature (list): A list containing which solutions at which steps should be used. Set to 1, for all times at
which you want to use the solution and to 0 at all other times
f_signature (list): Analogue to u_signature for the right hand side evaluations
checks (bool): Perform some checks on stability and convergence

Returns:
list: Coefficients for u
list: Coefficients for f
'''
n_u = np.sum(u_signature, dtype=int)
n_f = np.sum(f_signature, dtype=int)
n = n_u + n_f # number of available values
j = np.arange(n) # index variable
inv_fac = 1. / factorial(j) # compute inverse factorials once to use in Taylor expansions

# build a matrix containing the Taylor coefficients
A = np.zeros((n, n))

# fill the entries for u
for i in range(n_u):
A[:, i] = steps[u_signature > 0][i]**j * inv_fac

# fill the entries for f
for i in range(n_f):
A[1:, i + n_u] = steps[f_signature > 0][i]**j[:-1] * inv_fac[:-1]

# build a right hand side vector for solving the system
b = np.zeros(n)
b[0] = 1.

# solve the linear system
coeff = np.linalg.solve(A, b)

# distribute the coefficients
u_coeff = np.zeros_like(u_signature)
u_coeff[u_signature > 0] = coeff[0: n_u]
f_coeff = np.zeros_like(f_signature)
f_coeff[f_signature > 0] = coeff[n_u:]

if checks:
# check that the method is consistent
check_linear_difference_operator(u_coeff, f_coeff, steps)

# check if our method is zero stable
verify_root_condition(first_characteristic_polynomial(u_coeff))

# Check if the method is stable for h*lambda=-1
p = stability_polynomial(u_coeff, f_coeff, -1.)
strict_root_condition(p)

return u_coeff, f_coeff


def first_characteristic_polynomial(u_coeff, r=1):
'''
The first characteristic polynomial of a linear multistep method is equal to the coefficients of umultiplied with
powers of r.

Args:
u_coeff: The alpha coefficients for u of the LMM in order of descending time difference to the solution we want
r: The variable of the polynomial

Returns:
numpy.ndarray: List containing the polynomial in r. Set r=1 to get the coefficients.
'''
j = np.arange(len(u_coeff))
rho = np.zeros_like(u_coeff)
rho = -u_coeff * r**j
rho[-1] = r**len(u_coeff)
return rho[::-1] # reverse the order to go along with usual definitions


def second_characteristic_polynomial(f_coeff, r=1):
'''
The second characteristic polynomial of a linear multistep method is equal to the coefficients multiplied with
powers of r.

Args:
f_coeff: The alpha coefficients for f of the LMM in order of descending time difference to the solution we want
r: The variable of the polynomial

Returns:
numpy.ndarray: List containing the polynomial in r. Set r=1 to get the coefficients.
'''
j = np.arange(len(f_coeff))
sigma = np.zeros_like(f_coeff)
sigma = f_coeff * r**j
return sigma[::-1] # reverse the order to go along with usual definitions


def verify_root_condition(rho):
'''
For a linear multistep method to be convergent, we require that all roots of the first characteristic polynomial
are distinct and have modulus smaller or equal to one.
This root condition implies that the method is zero stable and Dahlquist's theorem states that a zero stable and
consistent method is convergent. If we can also show that the method is consistent, we have thus shown it is
convergent.

Args:
rho (numpy.ndarray): Coefficients of the first characteristic polynomial

Returns:
bool: Whether the root condition is satisfied.
'''
# compute the roots of the polynomial
roots = np.roots(rho)

# check the conditions
roots_distinct = len(np.unique(roots)) == len(roots)
# give some leeway because we introduce some numerical error when computing the roots
modulus_condition = all(abs(roots) <= 1. + 10. * np.finfo(float).eps)

# raise errors if we violate one of the conditions
assert roots_distinct, "Not all roots of the first characteristic polynomial of the LMM are distinct!"
assert modulus_condition, "Some of the roots of the first characteristic polynomial of the LMM have modulus larger \
one!"
return roots_distinct and modulus_condition


def check_linear_difference_operator(u_coeff, f_coeff, steps):
'''
Check if the linear multistep method is consistent by doing a Taylor expansion and testing if all terms cancel
except for the first, which should be one.

Args:
u_coeff (numpy.ndarray): Coefficients for u in the LMM
f_coeff (numpy.ndarray): Coefficients for f in the LMM
steps (numpy.ndarray): Steps from point of expansion

Returns:
None
'''
order = len(steps)
taylor_coeffs = np.zeros((len(u_coeff) + len(f_coeff), order))

# fill in the coefficients
for i in range(order):
# get expansions of u
if u_coeff[i] != 0:
taylor_coeffs[i, :] = u_coeff[i] * taylor_expansion(steps[i], order)

# get expansions of f
if f_coeff[i] != 0:
taylor_coeffs[order + i, 1:] = f_coeff[i] * taylor_expansion(steps[i], order - 1)

# check that all is well
L = np.sum(taylor_coeffs, axis=0)
want = np.zeros_like(L)
want[0] = 1.
assert all(np.isclose(L, want)), "Some derivatives do not cancel in the Taylor expansion!"

return None


def taylor_expansion(step, order):
'''
Get coefficients of a Taylor expansion.

Args:
step (float): Time difference from point around which we expand
order (int): The order up to which we want to expand

Returns:
numpy.ndarray: List containing the coefficients of the derivatives of u in the Taylor expansion
'''
j = np.arange(order)
return step**j / factorial(j)


def stability_polynomial(u_coeff, f_coeff, h_hat):
'''
Construct the stability polynomial for a value of h_hat = h * lambda.

Args:
u_coeff (numpy.ndarray): Coefficients for u in the LMM
f_coeff (numpy.ndarray): Coefficients for f in the LMM
h_hat (float) Parameter where you want to check stability

Returns:
numpy.ndarray: List containing the coefficients of the stability polynomial
'''
rho = first_characteristic_polynomial(u_coeff)
sigma = second_characteristic_polynomial(f_coeff)
return rho - h_hat * sigma


def strict_root_condition(p):
'''
Check whether the roots of the polynomial ae strictly smaller than one.

Args:
p (numpy.ndarray): Coefficients for the polynomial

Returns:
None
'''
roots = np.roots(p)
assert all(abs(roots) < 1), "Polynomial does not satisfy strict root condition!"
return None
Original file line number Diff line number Diff line change
Expand Up @@ -167,12 +167,13 @@ def solve_system(self, rhs, factor, u0, t):
me[:] = L.solve(rhs)
return me

def u_exact(self, t):
def u_exact(self, t, u_init=None, t_init=None):
"""
Routine to compute the exact solution at time t

Args:
t (float): current time
u_init, t_init: unused parameters for common interface reasons

Returns:
dtype_u: exact solution
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -169,12 +169,13 @@ def solve_system(self, rhs, factor, u0, t):
tol=self.params.lintol, maxiter=self.params.liniter)[0].reshape(self.params.nvars)
return me

def u_exact(self, t):
def u_exact(self, t, u_init=None, t_init=None):
"""
Routine to compute the exact solution at time t

Args:
t (float): current time
u_init, t_init: Dummy variables to provide unified interface

Returns:
dtype_u: exact solution
Expand Down
Loading