"""
Poisson equation: -\Delta u = f
"""
import numpy as np
from scipy import sparse
from typing import Dict, Any
from matrix_toolkit.pde.base import BasePDEEquation
from matrix_toolkit.pde.config import PDEConfig
[docs]
class PoissonEquation(BasePDEEquation):
"""
The Poisson equation:
.. math::
-\Delta u = f
In 1D:
.. math::
-\frac{d^2 u}{dx^2} = f
In 2D:
.. math::
-\left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \right) = f
In 3D:
.. math::
-\left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} + \frac{\partial^2 u}{\partial z^2} \right) = f
Examples:
>>> config = PDEConfig(dimension=2, mesh_size=64)
>>> eq = PoissonEquation(config)
>>> A = eq.generate_matrix()
"""
[docs]
def get_operator_name(self) -> str:
return "Laplacian"
[docs]
def get_default_coefficients(self) -> Dict[str, Any]:
return {
'diffusion': 1.0, # Coefficient of -Δ
}
[docs]
def randomize_coefficients(self, seed=None) -> Dict[str, Any]:
"""Randomize diffusion coefficient"""
coeffs = super().randomize_coefficients(seed)
# Ensure diffusion is actually randomized
if 'diffusion' in self.config.param_ranges:
if seed is not None:
np.random.seed(seed)
low, high = self.config.param_ranges['diffusion']
coeffs['diffusion'] = np.random.uniform(low, high)
return coeffs
[docs]
def generate_matrix(self) -> sparse.spmatrix:
"""Generate Poisson matrix for configured dimension"""
# Get coefficients
coeffs = self.config.coefficients.copy() if self.config.coefficients else self.get_default_coefficients()
# Apply randomization if needed
if self.config.random_params:
random_coeffs = self.randomize_coefficients(self.config.random_seed)
coeffs.update(random_coeffs)
diffusion = coeffs.get('diffusion', 1.0)
# Generate based on dimension
if self.config.dimension == 1:
return self._generate_1d(diffusion)
elif self.config.dimension == 2:
return self._generate_2d(diffusion)
elif self.config.dimension == 3:
return self._generate_3d(diffusion)
else:
raise ValueError(f"Unsupported dimension: {self.config.dimension}")
def _generate_1d(self, diffusion: float) -> sparse.spmatrix:
"""Generate 1D Poisson matrix"""
from matrix_toolkit.pde.dim1d.generators import generate_1d_laplacian
n = self.config.mesh_size[0]
h = self.config.get_mesh_spacing()[0]
A = generate_1d_laplacian(
n, h,
boundary=self.config.boundary_condition,
stencil_order=self.config.stencil_order
)
return diffusion * A
def _generate_2d(self, diffusion: float) -> sparse.spmatrix:
"""Generate 2D Poisson matrix"""
from matrix_toolkit.pde.dim2d.generators import generate_2d_laplacian
nx, ny = self.config.mesh_size
hx, hy = self.config.get_mesh_spacing()
A = generate_2d_laplacian(
nx, ny, hx, hy,
boundary=self.config.boundary_condition
)
return diffusion * A
def _generate_3d(self, diffusion: float) -> sparse.spmatrix:
"""Generate 3D Poisson matrix"""
from matrix_toolkit.pde.dim3d.generators import generate_3d_laplacian
nx, ny, nz = self.config.mesh_size
hx, hy, hz = self.config.get_mesh_spacing()
A = generate_3d_laplacian(
nx, ny, nz, hx, hy, hz,
boundary=self.config.boundary_condition
)
return diffusion * A
[docs]
def generate_rhs(self, source_function=None) -> np.ndarray:
"""
Generate right-hand side for Poisson equation
Args:
source_function: f(x) or f(x,y) or f(x,y,z)
Returns:
RHS vector
"""
if source_function is None:
# Default: f = 1
return np.ones(self.config.total_dofs())
# Create mesh grid
if self.config.dimension == 1:
x = np.linspace(
self.config.domain[0],
self.config.domain[1],
self.config.mesh_size[0]
)
return source_function(x)
elif self.config.dimension == 2:
nx, ny = self.config.mesh_size
x = np.linspace(self.config.domain[0], self.config.domain[1], nx)
y = np.linspace(self.config.domain[2], self.config.domain[3], ny)
X, Y = np.meshgrid(x, y, indexing='ij')
f = source_function(X, Y)
return f.flatten()
elif self.config.dimension == 3:
nx, ny, nz = self.config.mesh_size
x = np.linspace(self.config.domain[0], self.config.domain[1], nx)
y = np.linspace(self.config.domain[2], self.config.domain[3], ny)
z = np.linspace(self.config.domain[4], self.config.domain[5], nz)
X, Y, Z = np.meshgrid(x, y, z, indexing='ij')
f = source_function(X, Y, Z)
return f.flatten()