"""
3D PDE matrix generators
"""
import numpy as np
from scipy import sparse
from typing import Optional, Tuple
from matrix_toolkit.pde.config import BoundaryCondition
[docs]
def generate_3d_laplacian(
nx: int,
ny: int,
nz: int,
hx: float = None,
hy: float = None,
hz: float = None,
boundary: BoundaryCondition = BoundaryCondition.DIRICHLET
) -> sparse.spmatrix:
"""
Generate 3D Laplacian matrix: -(∂²u/∂x² + ∂²u/∂y² + ∂²u/∂z²)
Uses standard 7-point stencil.
Args:
nx, ny, nz: Grid points in each direction
hx, hy, hz: Grid spacings
boundary: Boundary condition
Returns:
Sparse CSR matrix of size (nx*ny*nz, nx*ny*nz)
Examples:
>>> A = generate_3d_laplacian(16, 16, 16)
>>> A.shape
(4096, 4096)
>>> A.nnz # 7-point stencil
28672
"""
if hx is None:
hx = 1.0 / (nx + 1)
if hy is None:
hy = 1.0 / (ny + 1)
if hz is None:
hz = 1.0 / (nz + 1)
# Create 1D Laplacians
from matrix_toolkit.pde.dim1d.generators import generate_1d_laplacian
Lx = generate_1d_laplacian(nx, hx, boundary)
Ly = generate_1d_laplacian(ny, hy, boundary)
Lz = generate_1d_laplacian(nz, hz, boundary)
# Identity matrices
Ix = sparse.eye(nx, format='csr')
Iy = sparse.eye(ny, format='csr')
Iz = sparse.eye(nz, format='csr')
# Kronecker sum: Lx ⊗ Iy ⊗ Iz + Ix ⊗ Ly ⊗ Iz + Ix ⊗ Iy ⊗ Lz
A = sparse.kron(sparse.kron(Lx, Iy), Iz, format='csr') + \
sparse.kron(sparse.kron(Ix, Ly), Iz, format='csr') + \
sparse.kron(sparse.kron(Ix, Iy), Lz, format='csr')
return A
[docs]
def generate_3d_convection_diffusion(
nx: int,
ny: int,
nz: int,
velocity: Tuple[float, float, float] = (1.0, 1.0, 1.0),
diffusion: float = 1.0,
hx: float = None,
hy: float = None,
hz: float = None
) -> sparse.spmatrix:
"""
Generate 3D convection-diffusion matrix
-ε∇²u + v·∇u
Args:
nx, ny, nz: Grid points
velocity: Convection velocity (vx, vy, vz)
diffusion: Diffusion coefficient
hx, hy, hz: Grid spacings
Returns:
Sparse CSR matrix
"""
if hx is None:
hx = 1.0 / (nx + 1)
if hy is None:
hy = 1.0 / (ny + 1)
if hz is None:
hz = 1.0 / (nz + 1)
vx, vy, vz = velocity
# Diffusion
A_diff = diffusion * generate_3d_laplacian(nx, ny, nz, hx, hy, hz)
# Convection derivatives
Dx = sparse.diags(
[-np.ones(nx-1), np.zeros(nx), np.ones(nx-1)],
[-1, 0, 1],
shape=(nx, nx)
) / (2 * hx)
Dy = sparse.diags(
[-np.ones(ny-1), np.zeros(ny), np.ones(ny-1)],
[-1, 0, 1],
shape=(ny, ny)
) / (2 * hy)
Dz = sparse.diags(
[-np.ones(nz-1), np.zeros(nz), np.ones(nz-1)],
[-1, 0, 1],
shape=(nz, nz)
) / (2 * hz)
Ix = sparse.eye(nx)
Iy = sparse.eye(ny)
Iz = sparse.eye(nz)
# Build convection operator
A_conv = vx * sparse.kron(sparse.kron(Dx, Iy), Iz, format='csr') + \
vy * sparse.kron(sparse.kron(Ix, Dy), Iz, format='csr') + \
vz * sparse.kron(sparse.kron(Ix, Iy), Dz, format='csr')
A = A_diff + A_conv
return A
[docs]
def generate_3d_heat(
nx: int,
ny: int,
nz: int,
dt: float,
diffusion: float = 1.0,
hx: float = None,
hy: float = None,
hz: float = None,
theta: float = 0.5
) -> sparse.spmatrix:
"""
Generate 3D heat equation matrix: ∂u/∂t = ε∇²u
Args:
nx, ny, nz: Grid points
dt: Time step
diffusion: Thermal diffusivity
hx, hy, hz: Grid spacings
theta: Implicitness parameter
Returns:
System matrix M
"""
if hx is None:
hx = 1.0 / (nx + 1)
if hy is None:
hy = 1.0 / (ny + 1)
if hz is None:
hz = 1.0 / (nz + 1)
# Laplacian
L = generate_3d_laplacian(nx, ny, nz, hx, hy, hz)
# Identity
n = nx * ny * nz
I = sparse.eye(n, format='csr')
# Theta-method
M = I - theta * dt * diffusion * L
return M
[docs]
def generate_3d_wave(
nx: int,
ny: int,
nz: int,
dt: float,
wave_speed: float = 1.0,
hx: float = None,
hy: float = None,
hz: float = None
) -> sparse.spmatrix:
"""
Generate 3D wave equation matrix: ∂²u/∂t² = c²∇²u
Args:
nx, ny, nz: Grid points
dt: Time step
wave_speed: Wave speed c
hx, hy, hz: Grid spacings
Returns:
System matrix M
"""
if hx is None:
hx = 1.0 / (nx + 1)
if hy is None:
hy = 1.0 / (ny + 1)
if hz is None:
hz = 1.0 / (nz + 1)
# Laplacian
L = generate_3d_laplacian(nx, ny, nz, hx, hy, hz)
# Identity
n = nx * ny * nz
I = sparse.eye(n, format='csr')
# Wave equation
M = I - (wave_speed * dt)**2 * L
return M
[docs]
def generate_3d_helmholtz(
nx: int,
ny: int,
nz: int,
wavenumber: float,
hx: float = None,
hy: float = None,
hz: float = None
) -> sparse.spmatrix:
"""
Generate 3D Helmholtz matrix: -∇²u - k²u
Args:
nx, ny, nz: Grid points
wavenumber: Wave number k
hx, hy, hz: Grid spacings
Returns:
Sparse CSR matrix
"""
# Laplacian
A = generate_3d_laplacian(nx, ny, nz, hx, hy, hz)
# Add -k² term
n = nx * ny * nz
A = A - (wavenumber**2) * sparse.eye(n, format='csr')
return A
[docs]
def generate_3d_biharmonic(
nx: int,
ny: int,
nz: int,
hx: float = None,
hy: float = None,
hz: float = None
) -> sparse.spmatrix:
"""
Generate 3D biharmonic matrix: ∇⁴u = ∇²∇²u
Args:
nx, ny, nz: Grid points
hx, hy, hz: Grid spacings
Returns:
Sparse CSR matrix
"""
# Compute Laplacian
L = generate_3d_laplacian(nx, ny, nz, hx, hy, hz)
# Biharmonic = Laplacian²
A = L @ L
return A