Unverified Commit 163f72ed authored by Robert Speck's avatar Robert Speck Committed by GitHub
Browse files

Merge pull request #128 from brownbaerchen/FD_stencils

Function to generate general finite difference stencils.
parents 3e7d21ed 56e41ee8
import logging
import numpy as np
from scipy.special import factorial
from pySDC.helpers.pysdc_helper import FrozenClass
......@@ -57,3 +59,42 @@ class ptype(object):
Abstract interface to apply mass matrix (only needed for FEM)
"""
raise NotImplementedError('ERROR: if you want a mass matrix, implement apply_mass_matrix(u)')
def get_finite_difference_stencil(derivative, order, type=None, steps=None):
"""
Derive general finite difference stencils from Taylor expansions
"""
if steps is not None:
n = len(steps)
elif type == 'center':
n = order + derivative - (derivative + 1) % 2 // 1
steps = np.arange(n) - (n) // 2
elif type == 'forward':
n = order + derivative
steps = np.arange(n)
elif type == 'backward':
n = order + derivative
steps = -np.arange(n)
else:
raise ValueError(f'Stencil must be of type "center", "forward" or "backward", not {type}. If you want something\
else, you can also give specific steps.')
# the index of the position around which we Taylor expand
zero_pos = np.argmin(abs(steps)) + 1
# make a matrix that contains the Taylor coefficients
A = np.zeros((n, n))
idx = np.arange(n)
inv_facs = 1. / factorial(idx)
for i in range(0, n):
A[i, :] = steps**idx[i] * inv_facs[i]
# make a right hand side vector that is zero everywhere except at the postition of the desired derivative
sol = np.zeros(n)
sol[derivative] = 1.
# solve the linear system for the finite difference coefficients
coeff = np.linalg.solve(A, sol)
return coeff, zero_pos, steps
......@@ -4,7 +4,7 @@ import scipy.sparse as sp
from scipy.sparse.linalg import cg
from pySDC.core.Errors import ParameterError, ProblemError
from pySDC.core.Problem import ptype
from pySDC.core.Problem import ptype, get_finite_difference_stencil
from pySDC.implementations.datatype_classes.mesh import mesh
......@@ -44,6 +44,8 @@ class heat2d_periodic(ptype):
raise ProblemError('need a square domain, got %s' % problem_params['nvars'])
if problem_params['nvars'][0] % 2 != 0:
raise ProblemError('the setup requires nvars = 2^p per dimension')
if 'order' not in problem_params:
problem_params['order'] = 2
# invoke super init, passing number of dofs, dtype_u and dtype_f
super(heat2d_periodic, self).__init__(init=(problem_params['nvars'], None, np.dtype('float64')),
......@@ -51,10 +53,10 @@ class heat2d_periodic(ptype):
# compute dx (equal in both dimensions) and get discretization matrix A
self.dx = 1.0 / self.params.nvars[0]
self.A = self.__get_A(self.params.nvars, self.params.nu, self.dx)
self.A = self.__get_A(self.params.nvars, self.params.nu, self.dx, self.params.order)
@staticmethod
def __get_A(N, nu, dx):
def __get_A(N, nu, dx, order):
"""
Helper function to assemble FD matrix A in sparse format
......@@ -67,16 +69,13 @@ class heat2d_periodic(ptype):
scipy.sparse.csc_matrix: matrix A in CSC format
"""
stencil = [1, -2, 1]
zero_pos = 2
stencil, zero_pos, _ = get_finite_difference_stencil(derivative=2, order=order, type='center')
dstencil = np.concatenate((stencil, np.delete(stencil, zero_pos - 1)))
offsets = np.concatenate(([N[0] - i - 1 for i in reversed(range(zero_pos - 1))],
[i - zero_pos + 1 for i in range(zero_pos - 1, len(stencil))]))
doffsets = np.concatenate((offsets, np.delete(offsets, zero_pos - 1) - N[0]))
A = sp.diags(dstencil, doffsets, shape=(N[0], N[0]), format='csc')
# stencil = [1, -2, 1]
# A = sp.diags(stencil, [-1, 0, 1], shape=(N[0], N[0]), format='csc')
A = sp.kron(A, sp.eye(N[0])) + sp.kron(sp.eye(N[1]), A)
A *= nu / (dx ** 2)
......
......@@ -23,15 +23,18 @@ def test_spatial_accuracy():
nvars_list = [(2 ** p, 2 ** p) for p in range(4, 12)]
# run accuracy test for all nvars
results = run_accuracy_check(nvars_list=nvars_list,problem_params=problem_params)
for order_stencil in [2, 4, 8]:
results = run_accuracy_check(nvars_list=nvars_list, problem_params=problem_params, order_stencil=order_stencil)
# compute order of accuracy
order = get_accuracy_order(results)
# compute order of accuracy
order = get_accuracy_order(results)
print(order_stencil, order)
assert (all(np.isclose(order, 2, rtol=0.005))), "ERROR: spatial order of accuracy is not as expected, got %s" %order
assert all(np.isclose(order, order_stencil, atol=5e-2)), f"ERROR: expected spatial order to be \
{order_stencil} but got {np.mean(order):.2f}"
def run_accuracy_check(nvars_list,problem_params):
def run_accuracy_check(nvars_list, problem_params, order_stencil):
"""
Routine to check the error of the Laplacian vs. its FD discretization
......@@ -49,6 +52,7 @@ def run_accuracy_check(nvars_list,problem_params):
# setup problem
problem_params['nvars'] = nvars
problem_params['order'] = order_stencil
prob = heat2d_periodic(problem_params=problem_params, dtype_u=mesh, dtype_f=mesh)
# create x values, use only inner points
......@@ -97,7 +101,8 @@ def get_accuracy_order(results):
id_prev = ID(nvars=nvars_list[i-1])
# compute order as log(prev_error/this_error)/log(this_nvars/old_nvars) <-- depends on the sorting of the list!
order.append(np.log(results[id_prev]/results[id])/np.log(nvars_list[i][0]/nvars_list[i-1][0]))
if results[id] > 1e-8 and results[id_prev] > 1e-8:
order.append(np.log(results[id_prev]/results[id])/np.log(nvars_list[i][0]/nvars_list[i-1][0]))
return order
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment