import numpy as np
import scipy.sparse as sp
from pySDC.implementations.problem_classes.boussinesq_helpers.buildFDMatrix import (
getMatrix,
getUpwindMatrix,
getBCLeft,
getBCRight,
)
#
#
#
[docs]
def get2DUpwindMatrix(N, dx, order):
Dx = getUpwindMatrix(N[0], dx, order)
return sp.kron(Dx, sp.eye(N[1]), format="csr")
#
#
#
[docs]
def get2DMesh(N, x_b, z_b, bc_hor, bc_ver):
assert np.size(N) == 2, 'N needs to be an array with two entries: N[0]=Nx and N[1]=Nz'
assert (
np.size(x_b) == 2
), 'x_b needs to be an array with two entries: x_b[0] = left boundary, x_b[1] = right boundary'
assert (
np.size(z_b) == 2
), 'z_b needs to be an array with two entries: z_b[0] = lower boundary, z_b[1] = upper boundary'
h = np.zeros(2)
x = None
z = None
if bc_hor[0] in ['periodic']:
assert bc_hor[1] in ['periodic'], 'Periodic boundary conditions must be prescribed at both boundaries'
x = np.linspace(x_b[0], x_b[1], N[0], endpoint=False)
h[0] = x[1] - x[0]
if bc_hor[0] in ['dirichlet', 'neumann']:
x = np.linspace(x_b[0], x_b[1], N[0] + 2, endpoint=True)
x = x[1 : N[0] + 1]
h[0] = x[1] - x[0]
if bc_ver[0] in ['periodic']:
assert bc_ver[1] in ['periodic'], 'Periodic boundary conditions must be prescribed at both boundaries'
z = np.linspace(z_b[0], z_b[1], N[1], endpoint=False)
h[1] = z[1] - z[0]
if bc_ver[0] in ['dirichlet', 'neumann']:
z = np.linspace(z_b[0], z_b[1], N[1] + 2, endpoint=True)
z = z[1 : N[1] + 1]
h[1] = z[1] - z[0]
xx, zz = np.meshgrid(x, z, indexing="ij")
return xx, zz, h
#
#
#
[docs]
def get2DMatrix(N, h, bc_hor, bc_ver, order):
assert np.size(N) == 2, 'N needs to be an array with two entries: N[0]=Nx and N[1]=Nz'
assert np.size(h) == 2, 'h needs to be an array with two entries: h[0]=dx and h[1]=dz'
Ax = getMatrix(N[0], h[0], bc_hor[0], bc_hor[1], order)
Az = getMatrix(N[1], h[1], bc_ver[0], bc_ver[1], order)
Dx = sp.kron(Ax, sp.eye(N[1]), format="csr")
Dz = sp.kron(sp.eye(N[0]), Az, format="csr")
return Dx, Dz
#
# NOTE: So far only constant dirichlet values can be prescribed, i.e. one fixed value for a whole segment
#
[docs]
def getBCHorizontal(value, N, dx, bc_hor):
assert (
np.size(value) == 2
), 'Value needs to be an array with two entries: value[0] for the left and value[1] for the right boundary'
assert np.size(N) == 2, 'N needs to be an array with two entries: N[0]=Nx and N[1]=Nz'
assert np.size(dx) == 1, 'dx must be a scalar'
assert (
np.size(bc_hor) == 2
), 'bc_hor must have two entries, bc_hor[0] specifying the BC at the left, bc_hor[1] at the right boundary'
bl = getBCLeft(value[0], N[0], dx, bc_hor[0])
bl = np.kron(bl, np.ones(N[1]))
br = getBCRight(value[1], N[0], dx, bc_hor[1])
br = np.kron(br, np.ones(N[1]))
return bl, br
[docs]
def getBCVertical(value, N, dz, bc_ver):
assert (
np.size(value) == 2
), 'Value needs to be an array with two entries: value[0] for the left and value[1] for the right boundary'
assert np.size(N) == 2, 'N needs to be an array with two entries: N[0]=Nx and N[1]=Nz'
assert np.size(dz) == 1, 'dx must be a scalar'
assert (
np.size(bc_ver) == 2
), 'bc_hor must have two entries, bc_hor[0] specifying the BC at the left, bc_hor[1] at the right boundary'
bd = getBCLeft(value[0], N[1], dz, bc_ver[0])
bd = np.kron(np.ones(N[0]), bd)
bu = getBCRight(value[1], N[1], dz, bc_ver[1])
bu = np.kron(np.ones(N[0]), bu)
return bd, bu