The Poisson equation connects a distributed source to a scalar potential whose gradient produces a field. In its most useful engineering form, $-\nabla\cdot(\kappa\nabla u)=f$ on a domain $\Omega$, it is not merely a formula containing a Laplacian: it is a boundary-value model whose coefficient, source, geometry, interfaces, gauge, and boundary conditions jointly define the answer. Electrostatic potential, steady temperature, gravitational potential, pressure correction, diffusion, and semiconductor band bending share this mathematical structure while assigning different units and signs to every symbol.

<svg viewBox="0 0 760 470" xmlns="http://www.w3.org/2000/svg" font-family="Arial,sans-serif"><rect x="0" y="0" width="760" height="470" rx="18" fill="#0d1117"/><text x="380" y="35" fill="#e6edf3" font-size="21" font-weight="700" text-anchor="middle">A Poisson problem is a complete source-to-response model</text><text x="380" y="62" fill="#8b949e" font-size="12" text-anchor="middle">Equation, material, domain, and boundary data are inseparable</text><rect x="28" y="125" width="150" height="190" rx="12" fill="#161b22" stroke="#f85149" stroke-width="2"/><text x="103" y="160" fill="#ff7b72" font-size="15" font-weight="700" text-anchor="middle">Source f</text><circle cx="78" cy="220" r="12" fill="#f85149"/><circle cx="128" cy="248" r="18" fill="#f85149" opacity=".65"/><text x="103" y="290" fill="#c9d1d9" font-size="11" text-anchor="middle">charge · heat · mass</text><path d="M178 220H220" stroke="#58a6ff" stroke-width="4"/><polygon points="220,220 207,212 207,228" fill="#58a6ff"/><rect x="220" y="100" width="320" height="240" rx="14" fill="#161b22" stroke="#58a6ff" stroke-width="2"/><text x="380" y="145" fill="#79c0ff" font-size="17" font-weight="700" text-anchor="middle">−∇·(κ∇u) = f in Ω</text><text x="380" y="195" fill="#e6edf3" font-size="12" text-anchor="middle">geometry Ω + coefficient κ(x)</text><text x="380" y="235" fill="#e6edf3" font-size="12" text-anchor="middle">Dirichlet / Neumann / Robin data</text><text x="380" y="275" fill="#e6edf3" font-size="12" text-anchor="middle">interfaces + reference condition</text><path d="M540 220H582" stroke="#3fb950" stroke-width="4"/><polygon points="582,220 569,212 569,228" fill="#3fb950"/><rect x="582" y="125" width="150" height="190" rx="12" fill="#161b22" stroke="#3fb950" stroke-width="2"/><text x="657" y="160" fill="#7ee787" font-size="15" font-weight="700" text-anchor="middle">Response u</text><path d="M610 270 Q657 150 704 270" fill="none" stroke="#3fb950" stroke-width="4"/><text x="657" y="290" fill="#c9d1d9" font-size="11" text-anchor="middle">potential and field</text><text x="380" y="410" fill="#c9d1d9" font-size="12" text-anchor="middle">A contour plot is meaningful only after every box in the model has been declared.</text></svg>

The divergence-form equation is the safest canonical statement. Writing $-\nabla\cdot(\kappa\nabla u)=f$ preserves conservation when $\kappa$ varies in space or jumps at an interface. Only for constant $\kappa$ may it be reduced without qualification to $-\kappa\nabla^2u=f$. Expanding a variable coefficient produces $-\kappa\nabla^2u-\nabla\kappa\cdot\nabla u=f$; silently dropping the second term changes the model. The coefficient may be scalar, a symmetric positive-definite tensor, or a nonlinear function of state. Uniform positive lower and upper bounds on it express ellipticity and prevent a nominal diffusion direction from disappearing.

Sign conventions must be fixed before sources are interpreted. Mathematicians often use $-\Delta u=f$ because the resulting Dirichlet operator is positive definite, while electrostatics is commonly written $\nabla\cdot(\epsilon\nabla\phi)=-\rho$. These are consistent under $u=\phi$, $\kappa=\epsilon$, and $f=\rho$. If electric field is $\mathbf E=-\nabla\phi$, Gauss's law $\nabla\cdot(\epsilon\mathbf E)=\rho$ recovers the same equation. A sign error reverses field direction, converts a potential minimum into a maximum, and can still leave a linear solver with a small residual, so sign tests belong in the physics verification rather than in solver diagnostics alone.

The Poisson equation is elliptic and normally poses a spatial boundary-value problem. Information is coupled across the whole connected domain: changing a boundary value or a source influences the solution everywhere, although Green-function influence may decay with distance. There is no initial condition or propagation speed in the static equation. A parabolic diffusion equation may relax toward a Poisson steady state, but pseudo-time marching is a numerical strategy, not the physical assertion that an electrostatic field slowly diffuses into equilibrium. This global character explains both the smoothness of solutions away from sources and the communication cost of scalable solvers.

Laplace's equation is the source-free special case rather than a different theory. Where $f=0$, the potential is harmonic and obeys the mean-value property: its value at a point equals an appropriate average over surrounding spheres or circles. Harmonic functions cannot have a strict interior maximum or minimum unless constant. Consequently, an observed interior electrostatic extremum in a charge-free, uniform-permittivity region signals a boundary effect, a coefficient interface, an unmodeled source, a numerical artifact, or a misread plotted quantity. Piecewise source-free regions can still have strongly curved global solutions because boundary and interface values transmit sources located elsewhere.

Dimensional analysis distinguishes physically similar-looking Poisson models. In SI electrostatics, $\epsilon$ has units of farads per meter, $\phi$ volts, and $\rho$ coulombs per cubic meter, so $\nabla\cdot(\epsilon\nabla\phi)$ has charge-density units. For steady heat conduction, $-\nabla\cdot(k\nabla T)=Q$ uses thermal conductivity and volumetric heat generation. Gravitational potential satisfies $\nabla^2\Phi=4\pi G\rho_m$. Pressure Poisson equations in incompressible flow inherit source units from velocity divergence and time step. Copying tolerances, scaling, or boundary interpretations across these applications without nondimensionalization is unsafe.

A complete model declares domain topology and coordinate system. Cartesian, cylindrical, and spherical Laplacians contain different geometric terms; an axisymmetric model has $r^{-1}\partial_r(r\partial_r u)$ and requires regularity on the axis. A two-dimensional cross-section may mean invariance per unit depth, axisymmetry, or a thin-sheet approximation, each assigning different source units. Holes and disconnected boundaries introduce additional harmonic degrees of freedom. Infinite domains need decay, radiation-like, periodic, or transformed boundary conditions rather than an arbitrary finite box whose edge quietly shapes the answer.

<svg viewBox="0 0 760 470" xmlns="http://www.w3.org/2000/svg" font-family="Arial,sans-serif"><rect x="0" y="0" width="760" height="470" rx="18" fill="#0d1117"/><text x="380" y="35" fill="#e6edf3" font-size="21" font-weight="700" text-anchor="middle">Boundary data choose one member of the solution family</text><text x="380" y="61" fill="#8b949e" font-size="12" text-anchor="middle">The same differential equation represents different physical systems</text><rect x="28" y="100" width="215" height="280" rx="14" fill="#161b22" stroke="#58a6ff" stroke-width="2"/><text x="135" y="135" fill="#79c0ff" font-size="15" font-weight="700" text-anchor="middle">Dirichlet</text><rect x="70" y="175" width="130" height="120" fill="#21262d" stroke="#58a6ff" stroke-width="5"/><text x="135" y="230" fill="#e6edf3" font-size="17" text-anchor="middle">u = g</text><text x="135" y="330" fill="#c9d1d9" font-size="11" text-anchor="middle">value imposed</text><rect x="273" y="100" width="215" height="280" rx="14" fill="#161b22" stroke="#3fb950" stroke-width="2"/><text x="380" y="135" fill="#7ee787" font-size="15" font-weight="700" text-anchor="middle">Neumann</text><rect x="315" y="175" width="130" height="120" fill="#21262d" stroke="#30363d" stroke-width="2"/><path d="M315 210H275M315 260H275M445 210H485M445 260H485" stroke="#3fb950" stroke-width="4"/><text x="380" y="230" fill="#e6edf3" font-size="15" text-anchor="middle">κ ∂ₙu = h</text><text x="380" y="330" fill="#c9d1d9" font-size="11" text-anchor="middle">normal flux imposed</text><rect x="518" y="100" width="215" height="280" rx="14" fill="#161b22" stroke="#d29922" stroke-width="2"/><text x="625" y="135" fill="#e3b341" font-size="15" font-weight="700" text-anchor="middle">Robin / mixed</text><rect x="560" y="175" width="130" height="120" fill="#21262d" stroke="#d29922" stroke-width="3"/><path d="M560 235H525" stroke="#d29922" stroke-width="4"/><text x="625" y="225" fill="#e6edf3" font-size="13" text-anchor="middle">αu + βκ∂ₙu</text><text x="625" y="250" fill="#e6edf3" font-size="13" text-anchor="middle">= r</text><text x="625" y="330" fill="#c9d1d9" font-size="11" text-anchor="middle">value and exchange coupled</text><text x="380" y="430" fill="#c9d1d9" font-size="11" text-anchor="middle">Pure Neumann data also require global compatibility and a gauge.</text></svg>

Dirichlet data prescribe the potential or primary field value. The condition $u=g$ on a boundary can represent a conductor held at known voltage, a surface held at known temperature, or a manufactured mathematical constraint. In finite elements it is an essential condition imposed on the trial space; in finite differences its known stencil contributions move to the right-hand side. A grounded electrostatic contact means a chosen reference voltage, not necessarily zero charge. Imposing measured potential along every accessible boundary may overconstrain a model if the measurement already includes contact drops not represented in the domain.

Neumann data prescribe normal flux and expose a compatibility condition. The condition $\kappa\partial_nu=h$ fixes outward conductive flux under the stated sign convention. Integrating the PDE gives $-\int_{\partial\Omega}\kappa\partial_nu\,dS=\int_\Omega f\,dV$, so source and boundary flux must balance for a pure-Neumann problem. If they do, adding any constant to $u$ gives another solution. One reference value, zero-mean constraint, or nullspace-aware solver fixes the gauge without inventing physics. If compatibility fails, a solver may stagnate, drift, or return a least-squares compromise rather than a genuine solution.

Robin conditions model exchange rather than an arbitrary algebraic mixture. A form such as $\alpha u+\beta\kappa\partial_nu=r$ can express convection to an ambient temperature, a surface capacitance, or a finite transfer impedance. The coefficients and outward-normal sign must match the physical balance. Taking limits recovers Dirichlet- or Neumann-like behavior, but extreme coefficient ratios can worsen conditioning. A boundary layer that was eliminated from the domain often leads to a Robin condition only after an asymptotic or circuit reduction; fitting it without that interpretation can hide frequency dependence and nonlocal effects.

Mixed boundaries partition the surface by physical role. A device may hold contact potentials on electrodes, impose symmetry flux on mirror planes, use insulating flux elsewhere, and approximate an open boundary on an exterior truncation. At corners where types meet, solution derivatives can be singular even when data are smooth. Those singularities reduce global convergence rates and motivate graded or adaptive meshes. Replacing every unspecified boundary by zero Neumann is not neutral: it asserts zero normal flux, and in electrostatics it can confine field lines that should leave the computational box.

Interface conditions follow from conservation and constitutive laws. In the absence of a singular sheet source, $u$ is normally continuous and normal flux $\kappa\partial_nu$ is continuous across an internal material boundary. A prescribed surface charge creates a flux jump, while an ideal dipole sheet may create a potential jump. Tangential derivatives follow from the trace where potential is continuous. Numerically smearing a sharp permittivity interface alters capacitance and field peaks; enforcing derivative continuity instead of displacement-flux continuity is wrong whenever coefficients differ.

The maximum principle provides a powerful qualitative audit. For $-\Delta u=f$ with $f\ge0$ under the chosen convention, curvature and extrema obey a constrained pattern; the exact inequality direction should be derived rather than memorized across sign conventions. With $f=0$, extrema occur at boundaries. Discrete matrices with suitable positive stencil structure inherit a discrete maximum principle. Overshoots in a nominally harmonic region can reveal a non-monotone high-order scheme, distorted finite elements, inconsistent boundary interpolation, or solver error. The principle does not apply unchanged to indefinite screened equations or nonlinear models.

Uniqueness comes from an energy identity. Subtract two solutions with the same Dirichlet data, multiply the homogeneous equation by their difference, integrate by parts, and obtain $\int_\Omega\kappa|\nabla w|^2dV=0$. Positive $\kappa$ forces $w$ constant, and a nonempty Dirichlet boundary fixes that constant to zero. The same reasoning shows why pure Neumann conditions retain exactly the constant nullspace on a connected domain. This proof is more than theory: it identifies the quadratic energy, predicts the algebraic matrix structure, and specifies which boundary choices should make conjugate gradients applicable.

The solution minimizes a potential-energy functional under Dirichlet constraints. For symmetric positive $\kappa$, $J[v]=\tfrac12\int\kappa|\nabla v|^2dV-\int fv\,dV$ has the Poisson solution as its minimizer over admissible fields. The first variation yields the weak equation. Convexity gives uniqueness once gauges are removed. Energy error is therefore a natural finite-element metric, and poor local gradients can matter even when pointwise potential appears smooth. If the physical model couples field energy to mobile charge, chemical free energy, or nonlinear polarization, the correct total functional may be nonlinear and the simple quadratic picture becomes only one block.

Weak solutions permit rough data and discontinuous material coefficients. Instead of demanding two classical derivatives everywhere, seek $u$ in an $H^1$ space so that $\int\kappa\nabla u\cdot\nabla v=\int fv$ plus boundary terms for all test functions $v$. This statement naturally encodes flux continuity and Neumann data. Lax–Milgram reasoning supplies existence and uniqueness when the bilinear form is coercive and data are bounded in the appropriate dual space. Point charges, reentrant corners, and abrupt interfaces reduce regularity, so expecting globally smooth second derivatives or textbook convergence rates can be mathematically unjustified.

<svg viewBox="0 0 760 470" xmlns="http://www.w3.org/2000/svg" font-family="Arial,sans-serif"><rect x="0" y="0" width="760" height="470" rx="18" fill="#0d1117"/><text x="380" y="35" fill="#e6edf3" font-size="21" font-weight="700" text-anchor="middle">A Green function maps each source point to its influence</text><text x="380" y="61" fill="#8b949e" font-size="12" text-anchor="middle">Geometry and boundary conditions are already built into G(x,ξ)</text><rect x="48" y="105" width="664" height="260" rx="16" fill="#161b22" stroke="#30363d" stroke-width="2"/><circle cx="225" cy="235" r="15" fill="#f85149"/><text x="225" y="205" fill="#ff7b72" font-size="13" text-anchor="middle">source ξ</text><circle cx="225" cy="235" r="45" fill="none" stroke="#d29922" stroke-width="2"/><circle cx="225" cy="235" r="85" fill="none" stroke="#d29922" stroke-width="2" opacity=".7"/><circle cx="225" cy="235" r="125" fill="none" stroke="#d29922" stroke-width="2" opacity=".4"/><circle cx="565" cy="190" r="8" fill="#58a6ff"/><text x="565" y="165" fill="#79c0ff" font-size="13" text-anchor="middle">observation x</text><path d="M240 230C350 195 430 190 557 190" fill="none" stroke="#3fb950" stroke-width="4" stroke-dasharray="8 5"/><text x="400" y="178" fill="#7ee787" font-size="13" text-anchor="middle">G(x,ξ)</text><path d="M48 105L712 365M48 365L712 105" stroke="#30363d" stroke-width="1" opacity=".4"/><text x="380" y="405" fill="#e6edf3" font-size="15" text-anchor="middle">u(x) = ∫Ω G(x,ξ) f(ξ) dξ + boundary contribution</text><text x="380" y="442" fill="#c9d1d9" font-size="11" text-anchor="middle">Free-space kernels alone do not satisfy a finite device boundary.</text></svg>

Green functions separate source superposition from boundary geometry. A Green function solves the operator equation for a unit point source with specified homogeneous boundary conditions. Then a distributed-source solution is an integral of $G(\mathbf x,\boldsymbol\xi)f(\boldsymbol\xi)$ plus the appropriate boundary contribution. Reciprocity $G(x,\xi)=G(\xi,x)$ follows for self-adjoint scalar problems. Changing the domain, coefficient, or boundary condition changes the Green function; using the free-space kernel inside a grounded enclosure without image or boundary corrections solves the wrong problem.

Fundamental solutions reveal dimension-dependent long-range behavior. For the Laplacian in three dimensions, the free-space kernel scales as $1/(4\pi r)$ up to the operator sign; in two dimensions it is logarithmic; in one dimension it is piecewise linear. These differences affect decay, neutrality requirements, and finite-domain sensitivity. A two-dimensional electrostatic point source actually represents an infinite line charge in a translationally invariant three-dimensional interpretation. Mesh-independent comparison therefore requires the correct dimensional source normalization rather than the same numerical delta value.

Distributional sources require integrated interpretation. A Dirac delta is not a large ordinary value at one grid node; it is defined by its action under integration. Depositing a point or particle charge onto a mesh should conserve total source and ideally preserve moments appropriate to the discretization. The exact continuum solution is singular, so pointwise error at the source does not converge in the usual sense. Regularization represents finite source size or numerical smoothing and changes self-energy and near-field peaks; its width must be reported and tested rather than hidden as a meshing detail.

The method of images is a boundary construction with limited geometries. Replacing a grounded plane or sphere by fictitious sources can reproduce the boundary condition in the physical region and yield exact fields. Image charges are not physical charges in the excluded conductor. The technique becomes cumbersome or unavailable for general shapes, multiple dielectric interfaces, and nonlinear media. It remains valuable as a verification case for numerical solvers because it supplies known forces, induced surface charge, and potential behavior near a boundary singularity.

Separation of variables turns simple boundaries into modal expansions. Rectangles, cylinders, and spheres admit eigenfunctions matched to coordinate surfaces, with Fourier, Bessel, or spherical-harmonic factors. A particular solution accounts for the source and a harmonic correction enforces boundary data. Series convergence can be slow near discontinuous boundary values and may exhibit Gibbs-like behavior in traces, while derivatives converge more delicately. Modal solutions provide reference values and scaling intuition, but complex device geometries usually require numerical discretization.

Multipole expansions compress distant source structure. Far from a localized distribution, total charge gives the monopole term, the first moment gives a dipole term, and higher moments decay more rapidly in three-dimensional free space. Neutral distributions have no monopole contribution. Boundaries and dielectric contrast modify this hierarchy through induced sources. Multipoles accelerate far-field evaluation and explain why local charge rearrangements can have weak distant influence, but truncation is controlled by separation ratio and fails when target and source regions overlap.

Boundary-integral methods move homogeneous-region work to surfaces. Green identities express a solution through boundary potential and normal derivative, reducing a three-dimensional homogeneous problem to a two-dimensional surface discretization. The resulting matrices are dense, singular quadrature requires care, and material junctions need consistent integral equations. Fast multipole and hierarchical methods can reduce cost. Boundary elements excel for open electrostatics with piecewise homogeneous media, while volumetric nonlinear charge or spatially varying coefficients often favor finite elements, finite volumes, or hybrid formulations.

<svg viewBox="0 0 760 470" xmlns="http://www.w3.org/2000/svg" font-family="Arial,sans-serif"><rect x="0" y="0" width="760" height="470" rx="18" fill="#0d1117"/><text x="380" y="35" fill="#e6edf3" font-size="21" font-weight="700" text-anchor="middle">Finite differences turn curvature into a sparse stencil</text><text x="380" y="61" fill="#8b949e" font-size="12" text-anchor="middle">Five local couplings assemble into one global elliptic system</text><g stroke="#30363d" stroke-width="1"><path d="M70 100V370M130 100V370M190 100V370M250 100V370M310 100V370M40 130H340M40 190H340M40 250H340M40 310H340"/></g><g fill="#58a6ff"><circle cx="190" cy="250" r="13"/><circle cx="130" cy="250" r="10"/><circle cx="250" cy="250" r="10"/><circle cx="190" cy="190" r="10"/><circle cx="190" cy="310" r="10"/></g><text x="190" y="255" fill="#0d1117" font-size="12" font-weight="700" text-anchor="middle">4</text><text x="130" y="254" fill="#0d1117" font-size="11" text-anchor="middle">−1</text><text x="250" y="254" fill="#0d1117" font-size="11" text-anchor="middle">−1</text><text x="190" y="194" fill="#0d1117" font-size="11" text-anchor="middle">−1</text><text x="190" y="314" fill="#0d1117" font-size="11" text-anchor="middle">−1</text><path d="M355 235H405" stroke="#3fb950" stroke-width="4"/><polygon points="405,235 392,227 392,243" fill="#3fb950"/><rect x="420" y="100" width="285" height="270" rx="14" fill="#161b22" stroke="#3fb950" stroke-width="2"/><g fill="#3fb950"><rect x="450" y="130" width="20" height="20"/><rect x="475" y="155" width="20" height="20"/><rect x="500" y="180" width="20" height="20"/><rect x="525" y="205" width="20" height="20"/><rect x="550" y="230" width="20" height="20"/><rect x="575" y="255" width="20" height="20"/><rect x="600" y="280" width="20" height="20"/><rect x="625" y="305" width="20" height="20"/></g><g fill="#58a6ff" opacity=".75"><rect x="475" y="130" width="20" height="20"/><rect x="450" y="155" width="20" height="20"/><rect x="500" y="155" width="20" height="20"/><rect x="475" y="180" width="20" height="20"/><rect x="525" y="180" width="20" height="20"/><rect x="500" y="205" width="20" height="20"/><rect x="550" y="205" width="20" height="20"/></g><text x="562" y="345" fill="#c9d1d9" font-size="12" text-anchor="middle">sparse matrix A</text><text x="380" y="420" fill="#e6edf3" font-size="14" text-anchor="middle">(4uᵢⱼ − uᵢ₋₁ⱼ − uᵢ₊₁ⱼ − uᵢⱼ₋₁ − uᵢⱼ₊₁)/h² = fᵢⱼ</text></svg>

The centered finite-difference stencil is simple only on a uniform Cartesian grid. In two dimensions, the five-point approximation to $-\Delta u$ is second-order accurate for a sufficiently smooth solution, and three dimensions gives a seven-point stencil. Taylor expansion establishes local truncation error, but global error also depends on boundaries, stability, source regularity, and geometry representation. A fine rectangular grid can be excellent for boxes and periodic cells. On curved domains, stair-step boundaries may dominate error even though the interior stencil remains formally second order.

Variable coefficients belong on stencil faces through fluxes. Discretize $-\nabla\cdot(\kappa\nabla u)$ by first estimating face flux and then differencing its divergence. Harmonic averaging is often appropriate for normal transport through layered media because it respects series resistance and flux continuity; arithmetic averaging can overpredict flux across a large contrast. Tensor coefficients create cross couplings and require schemes that preserve symmetry, consistency, and preferably monotonicity. Directly multiplying a constant-coefficient Laplacian row by nodal $\kappa$ generally fails to conserve interface flux.

Finite volumes make local conservation explicit. Integrating the PDE over each control volume gives a balance between outward face flux and integrated source. Neighboring cells share equal and opposite numerical flux, so global conservation follows by cancellation. Cell-centered methods handle material inventories naturally and are common in transport codes. Nonorthogonal meshes require gradient reconstruction and correction terms; aggressive corrections can lose monotonicity. Conservation does not by itself ensure accurate potential or gradient, so consistency and mesh-convergence tests remain necessary.

Finite elements begin from the weak form and accommodate complex geometry. Choose a mesh, an $H^1$-conforming approximation space, and basis functions; assemble stiffness entries $A_{ij}=\int\kappa\nabla N_i\cdot\nabla N_j$ and load entries from sources and natural boundaries. Piecewise linear elements on triangles or tetrahedra are robust, while higher order can converge rapidly for smooth geometry and fields. Curved boundaries should be represented at compatible order. MFEM's canonical Poisson examples make this sequence explicit: mesh, finite-element space, bilinear and linear forms, essential boundary elimination, solve, and error evaluation.

Essential and natural boundary conditions enter finite elements differently. Dirichlet data restrict trial degrees of freedom and require a consistent lifting for nonzero values. Neumann flux appears from integration by parts as a boundary integral and is therefore natural. Robin data contribute both matrix and load terms. Forgetting a natural boundary term implicitly imposes zero flux. Strongly setting nodal values on a curved or high-order boundary can reduce accuracy if the geometry mapping is inconsistent, while penalty or Nitsche enforcement introduces parameters and stability conditions that must be verified.

Mesh refinement must target the quantity of interest as well as the residual. Uniform $h$ refinement reduces cell size everywhere; $p$ refinement raises polynomial order; $hp$ strategies combine them. Residual-based estimators locate unresolved source, interface, corner, and boundary effects. Goal-oriented estimators emphasize an output such as capacitance or peak field rather than global energy error. Refining on a visually steep contour alone can waste elements, and refining around a mathematical point singularity may never make the pointwise peak converge. Report degrees of freedom, mesh family, estimator, and observed rate.

Discretization produces a sparse linear algebra problem whose structure carries physics. With Dirichlet anchoring and positive scalar coefficient, the stiffness matrix is normally symmetric positive definite. Pure Neumann data give a symmetric positive-semidefinite matrix with the constant vector in its nullspace. Nonsymmetric boundary treatments, advection-like couplings, multiphysics linearization, or poor elimination can change this classification. Solver choice should follow the actual assembled operator rather than the PDE's informal name. A matrix diagnostic can catch symmetry loss, null rows, coefficient sign errors, and disconnected unconstrained components before iteration begins.

Conditioning deteriorates as the mesh resolves smaller scales. For a basic Laplacian discretization, the condition number grows roughly like $h^{-2}$, and material contrast or stretched elements can make it worse. A small residual does not imply small solution error without considering conditioning and norm. Scaling unknowns and equations improves numerical representation but does not replace preconditioning. Double precision may be insufficient for extreme contrasts or nearly floating subdomains. Iteration histories should include the true residual and a physically meaningful stopping criterion, not only an implementation-specific preconditioned norm.

Conjugate gradients require a symmetric positive-definite operator. Applied to the anchored scalar Poisson matrix, CG builds energy-optimal approximations using short recurrences. SciPy's current sparse-linear-algebra documentation explicitly describes cg for symmetric positive-definite systems. It is inappropriate for an unprojected singular pure-Neumann matrix or a nonsymmetric assembly. MINRES can address symmetric indefinite or compatible singular cases under suitable handling; GMRES handles nonsymmetry at higher storage cost. Sparse direct solvers are valuable for modest problems and repeated right-hand sides but can suffer fill-in in three dimensions.

Preconditioning determines whether an iterative method scales. Jacobi rescales by the diagonal; incomplete factorization approximates elimination; algebraic multigrid builds a hierarchy from matrix connectivity; geometric multigrid uses known meshes. A useful preconditioner clusters eigenvalues or reduces error components cheaply, not necessarily approximates every entry. Strong coefficient jumps, anisotropy, thin layers, and mixed boundary conditions can defeat default coarsening or smoothers. Reusing a setup across nearby nonlinear iterations may save time, but changes in depletion, active regions, or coefficients can require rebuilding it.

<svg viewBox="0 0 760 470" xmlns="http://www.w3.org/2000/svg" font-family="Arial,sans-serif"><rect x="0" y="0" width="760" height="470" rx="18" fill="#0d1117"/><text x="380" y="35" fill="#e6edf3" font-size="21" font-weight="700" text-anchor="middle">Multigrid attacks error at the scale where it is inexpensive</text><text x="380" y="61" fill="#8b949e" font-size="12" text-anchor="middle">Smoothing and coarse-grid correction complement each other</text><rect x="40" y="100" width="205" height="250" rx="12" fill="#161b22" stroke="#58a6ff" stroke-width="2"/><text x="142" y="130" fill="#79c0ff" font-size="14" font-weight="700" text-anchor="middle">Fine grid</text><path d="M60 240Q75 150 90 240T120 240T150 240T180 240T225 240" fill="none" stroke="#f85149" stroke-width="3"/><text x="142" y="300" fill="#c9d1d9" font-size="11" text-anchor="middle">smooth high-frequency error</text><path d="M245 225H280" stroke="#3fb950" stroke-width="4"/><polygon points="280,225 267,217 267,233" fill="#3fb950"/><rect x="280" y="125" width="200" height="200" rx="12" fill="#161b22" stroke="#3fb950" stroke-width="2"/><text x="380" y="155" fill="#7ee787" font-size="14" font-weight="700" text-anchor="middle">Coarse grid</text><path d="M305 250Q340 165 375 250T455 250" fill="none" stroke="#d29922" stroke-width="4"/><text x="380" y="290" fill="#c9d1d9" font-size="11" text-anchor="middle">solve low-frequency error</text><path d="M480 225H515" stroke="#3fb950" stroke-width="4"/><polygon points="515,225 502,217 502,233" fill="#3fb950"/><rect x="515" y="100" width="205" height="250" rx="12" fill="#161b22" stroke="#a371f7" stroke-width="2"/><text x="617" y="130" fill="#d2a8ff" font-size="14" font-weight="700" text-anchor="middle">Correct and smooth</text><path d="M540 240Q585 215 630 240T695 240" fill="none" stroke="#a371f7" stroke-width="3"/><text x="617" y="300" fill="#c9d1d9" font-size="11" text-anchor="middle">mesh-independent convergence</text><text x="380" y="410" fill="#e6edf3" font-size="14" text-anchor="middle">restrict residual → coarse solve → prolong correction → post-smooth</text><text x="380" y="442" fill="#c9d1d9" font-size="11" text-anchor="middle">Coefficient-aware transfers are essential across strong material interfaces.</text></svg>

Multigrid resolves the frequency complementarity of elliptic error. Jacobi or Gauss–Seidel relaxation quickly damps error that oscillates from node to node but barely changes smooth error. On a coarser grid, that smooth error appears higher frequency and is cheap to correct. A V-cycle restricts residual, solves or relaxes on coarse levels, prolongs correction, and post-smooths. With appropriate transfer, smoothing, and coarse spaces, work can approach linear complexity in unknown count. Poor treatment of anisotropy or disconnected high-conductivity regions destroys this ideal behavior.

Fast transform solvers exploit separability rather than general sparsity. On rectangles with constant coefficients and compatible boundary conditions, discrete sine, cosine, or Fourier transforms diagonalize coordinate operators and produce very fast Poisson solves. Periodic problems require the zero Fourier mode to satisfy neutrality and be assigned a gauge. Embedded objects, irregular coefficients, and local refinement break direct separability, though FFT solvers may still serve as preconditioners. Padding and periodic images can contaminate open-domain electrostatics unless the Green kernel and cell size are designed for isolation.

Domain decomposition distributes work but creates interface obligations. Schwarz methods solve overlapping or nonoverlapping subdomains and exchange boundary information. Krylov methods with additive Schwarz preconditioners can scale across processors when a global coarse space communicates long-wavelength error. Without that coarse level, iteration count grows with subdomain count because the elliptic field is global. Load balance should include adaptive mesh density and nonlinear material work, not only cell count. Reproducibility can change slightly with reduction order, so tolerances should exceed floating-point communication noise.

Adaptive mesh refinement must preserve conservative transfer and solver hierarchy. Refining cells near junctions, corners, and localized charge reduces degrees of freedom relative to a uniformly fine mesh. Hanging-node constraints, coarse–fine flux consistency, and source projection need explicit treatment. When particles or dopants deposit charge, transferring between levels must conserve total charge. An adaptive sequence should demonstrate that the target output stabilizes and that refinement indicators shrink; a beautiful locally dense mesh is not evidence of accuracy by itself.

GPU acceleration rewards regular arithmetic but does not remove global coupling. Matrix-free stencil and high-order finite-element kernels can achieve high bandwidth and avoid storing sparse matrices. Krylov dot products require global reductions, triangular incomplete-factor solves offer limited parallelism, and coarse multigrid levels may underutilize a device. Mixed precision can accelerate smoothers or preconditioners while the outer residual is checked in higher precision. Performance reports need end-to-end setup, transfers, nonlinear iterations, and energy use, not kernel throughput alone.

Nondimensionalization exposes controlling ratios and improves numerical scale. Choose characteristic length $L$, potential $U$, coefficient $K$, and source $F$ so the normalized equation has order-one variables. Electrostatic semiconductor scaling may use thermal voltage $V_T=k_BT/q$ and Debye length, revealing stiffness when device and screening lengths differ greatly. Scaling does not change the physical solution when transformed back, but it makes tolerances comparable, protects exponential carrier laws from overflow, and clarifies which terms can be asymptotically neglected.

A residual is necessary but not sufficient evidence of a correct solution. The algebraic residual $r=b-Au_h$ measures satisfaction of the discrete equations. It does not measure truncation error, geometry error, incorrect coefficients, bad source units, or wrong boundary data. A solver can converge exactly to the wrong discretized model. Compare residual reduction with discretization estimates, global balance, mesh changes, and analytical limits. For nonlinear problems, distinguish the inner linear residual from the outer nonlinear residual and from the original dimensional PDE imbalance.

<svg viewBox="0 0 760 470" xmlns="http://www.w3.org/2000/svg" font-family="Arial,sans-serif"><rect x="0" y="0" width="760" height="470" rx="18" fill="#0d1117"/><text x="380" y="35" fill="#e6edf3" font-size="21" font-weight="700" text-anchor="middle">Semiconductor electrostatics closes a nonlinear charge loop</text><text x="380" y="61" fill="#8b949e" font-size="12" text-anchor="middle">Potential changes carriers, and carriers change potential</text><rect x="55" y="130" width="185" height="190" rx="14" fill="#161b22" stroke="#58a6ff" stroke-width="2"/><text x="147" y="165" fill="#79c0ff" font-size="14" font-weight="700" text-anchor="middle">Poisson solve</text><text x="147" y="210" fill="#e6edf3" font-size="13" text-anchor="middle">∇·(ε∇φ) = −ρ</text><text x="147" y="260" fill="#c9d1d9" font-size="11" text-anchor="middle">potential and field</text><path d="M240 225H290" stroke="#3fb950" stroke-width="4"/><polygon points="290,225 277,217 277,233" fill="#3fb950"/><rect x="290" y="130" width="180" height="190" rx="14" fill="#161b22" stroke="#3fb950" stroke-width="2"/><text x="380" y="165" fill="#7ee787" font-size="14" font-weight="700" text-anchor="middle">Carrier model</text><text x="380" y="210" fill="#e6edf3" font-size="12" text-anchor="middle">n(φ), p(φ), traps</text><text x="380" y="260" fill="#c9d1d9" font-size="11" text-anchor="middle">statistics or transport</text><path d="M470 225H520" stroke="#d29922" stroke-width="4"/><polygon points="520,225 507,217 507,233" fill="#d29922"/><rect x="520" y="130" width="185" height="190" rx="14" fill="#161b22" stroke="#d29922" stroke-width="2"/><text x="612" y="165" fill="#e3b341" font-size="14" font-weight="700" text-anchor="middle">Charge update</text><text x="612" y="210" fill="#e6edf3" font-size="12" text-anchor="middle">q(p−n+Nᴅ⁺−Nₐ⁻)</text><text x="612" y="260" fill="#c9d1d9" font-size="11" text-anchor="middle">mix and test convergence</text><path d="M612 320V375H147V320" fill="none" stroke="#a371f7" stroke-width="4"/><polygon points="147,320 139,333 155,333" fill="#a371f7"/><text x="380" y="400" fill="#d2a8ff" font-size="12" text-anchor="middle">Gummel iteration or coupled Newton method</text><text x="380" y="442" fill="#c9d1d9" font-size="11" text-anchor="middle">Converge potential, charge, current, and conserved terminal quantities together.</text></svg>

Semiconductor Poisson charge is state dependent. A common convention is $\rho=q(p-n+N_D^+-N_A^-)+\rho_{trap}+\rho_{fixed}$, with ionization, carrier statistics, traps, polarization, and fixed interface charge declared separately. Electron charge contributes negatively even though the elementary charge $q$ is positive. Doping is not always fully ionized, especially at low temperature or high degeneracy. Substituting a net-doping profile for total charge ignores mobile screening and is valid only in approximations such as selected depletion regions.

Band energies and electrostatic potential must share one energy reference. Electron potential energy changes as $-q\phi$, so raising electrostatic potential lowers electron band-edge energy under the usual convention. Work functions, electron affinity, band offsets, Fermi levels, and applied terminal voltages must be aligned consistently. A gauge shift in $\phi$ accompanied by the corresponding energy-reference shift changes no observable. Mixing electron-volts and volts without the factor $q$, or mixing electrostatic and electron potential signs, can produce plausible inverted band diagrams.

The depletion approximation is a controlled piecewise-charge model. In an abrupt pn junction, assume mobile carriers are negligible within depletion widths, leaving approximately $-qN_A$ on the p side and $+qN_D$ on the n side. Integrating Poisson twice with field continuity, charge neutrality $N_Ax_p=N_Dx_n$, and the built-in voltage yields the familiar square-root depletion width. MIT material examples use this structure for GaAs junction electrostatics. The approximation weakens near depletion edges, under high injection, in graded doping, nanoscale confinement, and when traps or incomplete ionization matter.

A MOS capacitor couples oxide Laplace behavior to semiconductor Poisson behavior. Ideal oxide has negligible mobile volume charge, so potential is linear in a one-dimensional uniform oxide, while semiconductor charge bends the bands nonlinearly. Gate work-function difference, oxide thickness and permittivity, fixed oxide charge, interface traps, substrate doping, and temperature set the voltage partition. Accumulation, depletion, and inversion are regimes of the same boundary-value problem. Treating the gate voltage as semiconductor surface potential discards oxide drop and charge and gives incorrect threshold and capacitance.

Heterojunctions demand displacement continuity and band-offset bookkeeping. Permittivity may jump while normal electric displacement remains continuous unless sheet charge is present. Electrostatic potential is generally continuous across a conventional interface, while conduction and valence band edges have material offsets beyond electrostatic bending. Polarization sheets in III-nitrides deliberately create displacement jumps and high-density channels. Averaging permittivity or doping across an interface can blur sheet physics, shift confined charge, and corrupt capacitance. Mesh faces should align with sharp interfaces when possible.

Drift–diffusion makes Poisson one equation in a coupled conservation system. Electron and hole continuity equations determine nonequilibrium carrier densities and currents, while Poisson supplies the electric field. Gummel iteration alternates subproblems with damping; Newton methods linearize the coupled residual for faster local convergence but require accurate Jacobians and robust globalization. Convergence of potential alone is insufficient if terminal currents or continuity residuals still change. High fields may also require mobility, generation–recombination, impact ionization, or energy-transport models whose validity bounds should be stated.

Poisson–Boltzmann equations are nonlinear through exponential populations. In a nondegenerate semiconductor or electrolyte, mobile species follow Boltzmann factors in potential, producing a nonlinear source. Linearization gives a screened Poisson or Debye–Hückel equation only when dimensionless potential is small. Large surface potentials, multivalent ions, steric crowding, degeneracy, and correlations violate that limit. Exponentials can overflow numerically; nondimensional variables, bounded line searches, continuation in bias, and stable carrier evaluations are practical necessities rather than cosmetic implementation choices.

Poisson–Schrödinger coupling turns quantum density into electrostatic source. Solve Schrödinger's equation in a confinement potential derived partly from $\phi$, occupy its states according to reservoirs and statistics, construct carrier density, and return that density to Poisson. This loop captures subbands and wavefunction penetration that classical local carrier laws miss. Exchange-correlation, effective masses, valley degeneracy, open boundaries, and nonequilibrium occupation determine the model. Mixing or Newton-like acceleration is usually required because an electrostatic update shifts the very levels that set charge.

Plasma sheaths use Poisson with kinetic or fluid charge closure. Electron and ion densities differ near material boundaries, creating a sheath electric field. Boltzmann electrons plus an ion model can yield a nonlinear sheath equation, while kinetic simulations deposit particle charge and solve Poisson repeatedly. Debye length dictates spatial resolution, and wall potential or current balance supplies boundary physics. Quasineutral bulk models deliberately omit small charge separation and cannot resolve a sheath without matching. This page's general framework complements, rather than replaces, dedicated plasma-chamber modeling.

Electrostatic capacitance is a derivative of charge with respect to voltage. Solve a sequence of linear Laplace or Poisson problems with conductor boundary potentials, integrate normal displacement to obtain terminal charge, and form a capacitance matrix or differential capacitance. The matrix should respect charge conservation and reciprocity under the model assumptions. In nonlinear semiconductors capacitance depends on bias, frequency, carrier response time, and trap kinetics; a static derivative is not automatically the measured high-frequency C–V curve. Numerical differencing step must exceed solver noise yet remain locally linear.

<svg viewBox="0 0 760 470" xmlns="http://www.w3.org/2000/svg" font-family="Arial,sans-serif"><rect x="0" y="0" width="760" height="470" rx="18" fill="#0d1117"/><text x="380" y="35" fill="#e6edf3" font-size="21" font-weight="700" text-anchor="middle">Verification triangulates equation, discretization, and physics</text><text x="380" y="61" fill="#8b949e" font-size="12" text-anchor="middle">No single residual or contour plot closes the evidence chain</text><polygon points="380,95 675,345 85,345" fill="#161b22" stroke="#58a6ff" stroke-width="3"/><circle cx="380" cy="120" r="62" fill="#21262d" stroke="#3fb950" stroke-width="3"/><text x="380" y="112" fill="#7ee787" font-size="14" font-weight="700" text-anchor="middle">Manufactured</text><text x="380" y="135" fill="#c9d1d9" font-size="11" text-anchor="middle">solution and rate</text><circle cx="140" cy="320" r="62" fill="#21262d" stroke="#d29922" stroke-width="3"/><text x="140" y="312" fill="#e3b341" font-size="14" font-weight="700" text-anchor="middle">Conservation</text><text x="140" y="335" fill="#c9d1d9" font-size="11" text-anchor="middle">source = flux</text><circle cx="620" cy="320" r="62" fill="#21262d" stroke="#a371f7" stroke-width="3"/><text x="620" y="312" fill="#d2a8ff" font-size="14" font-weight="700" text-anchor="middle">Benchmark</text><text x="620" y="335" fill="#c9d1d9" font-size="11" text-anchor="middle">analytic or measured</text><text x="380" y="275" fill="#e6edf3" font-size="15" font-weight="700" text-anchor="middle">credible field</text><text x="380" y="410" fill="#c9d1d9" font-size="12" text-anchor="middle">Also test gauge, symmetry, maximum principle, interfaces, and mesh independence.</text><text x="380" y="442" fill="#8b949e" font-size="11" text-anchor="middle">Validation asks whether the chosen Poisson model represents the experiment.</text></svg>

Manufactured solutions verify implementation without needing a natural exact case. Choose a smooth $u_{exact}$ compatible with the geometry, apply the differential operator to generate $f$, and derive boundary data from the same field. Solve on a mesh sequence and measure error in $L^2$, gradient, and outputs. The observed rate should match element or stencil theory until roundoff or solver tolerance dominates. Manufactured tests should exercise variable coefficients, each boundary type, curved geometry, and interfaces; a single constant-coefficient Dirichlet box leaves major code paths untested.

Global flux balance is the discrete form of Gauss's law. Integrate the numerical source and compare it with signed boundary flux, including sheet sources and electrode charge. Finite volumes may satisfy this locally by construction; finite elements satisfy weak balances whose evaluation requires consistent numerical flux. A small mismatch can be normalized by total absolute source or a relevant terminal quantity. Exact global balance can coexist with locally wrong fields, but failure of balance immediately identifies sign, boundary-normal, source-deposition, or nonlinear-convergence defects.

Mesh-convergence studies must hold the physical model fixed. Refine geometry, coefficients, sources, and boundaries consistently while tightening the algebraic tolerance enough that discretization dominates. Compare at least three credible resolutions and estimate observed order in smooth regimes. Peak field at a sharp metal corner may diverge with refinement, so use an averaged field, energy, force, or rounded physical geometry instead. If a model parameter such as interface width changes with cell size, the study measures a changing model rather than numerical convergence.

Analytical limits expose errors that residual checks cannot see. Recover a linear potential for one-dimensional source-free uniform material, a parabola for constant source, radial $1/r$ behavior outside a spherical source, charge neutrality for pure periodic domains, and the depletion-width scaling of an abrupt junction. Symmetry planes should show zero normal flux. Superposition should hold for a linear model. Translating every prescribed potential by a constant should leave fields unchanged when the gauge permits it. These tests are inexpensive and directly tied to governing structure.

Validation compares observable predictions to experiments through a measurement model. Electrostatic potential inside a device is rarely measured directly. C–V, Kelvin probe, electron holography, junction depletion, terminal charge, force, temperature, or current each filters the field through instrument response and uncertain geometry. Calibrating unknown fixed charge against one curve and validating against the same curve is not independent evidence. Report uncertainty in permittivity, doping, dimensions, contact work functions, traps, and temperature, then test predictions across conditions not used for fitting.

Sensitivity is often concentrated at boundaries and interfaces. A small work-function shift, oxide thickness error, surface-charge density, corner radius, or outer-domain location can dominate a well-converged interior discretization. Adjoint methods compute derivatives of a scalar output with respect to many parameters at cost comparable to a few solves. Sensitivities are local to the assumed model and do not quantify structural uncertainty from missing physics. Nondifferentiable mesh changes and nonlinear branch switching need special care.

Inverse Poisson problems are generally more fragile than forward solves. Inferring source $f$ from noisy potential requires differentiation, which amplifies high-frequency noise. Inferring permittivity, boundary charge, or geometry can be nonunique because different causes produce similar fields. Regularization, prior information, sensor models, and identifiability analysis are essential. A smooth reconstructed charge map may reflect the regularizer more than the data. Forward-solver accuracy should be substantially tighter than data discrepancy so numerical error is not misinterpreted as inferred structure.

Reproducibility requires recording the mathematical problem, not only software settings. Preserve the exact domain and units, coordinate interpretation, coefficient fields, source normalization, boundary partition and outward signs, interface laws, gauge, mesh, discretization order, linear and nonlinear tolerances, solver and preconditioner versions, and postprocessing definitions. A screenshot and a mesh count cannot reconstruct a calculation. Hashing input fields and exporting integrated balance diagnostics makes later comparisons more reliable than relying on mutable project files.

A practical workflow begins with conservation and ends with an observable. State what is sourced, what flux transports it, and what potential drives that flux. Fix units and signs, sketch material and boundary regions, identify gauges and compatibility, then select a discretization whose conservation and geometry properties match the problem. Classify the assembled operator before choosing a solver. Verify against a manufactured case, global balance, analytical limit, and mesh sequence. Only then compare an output defined through the actual measurement or design decision.

Define physical source f, coefficient κ, domain Ω, units, and sign convention
  -> Partition boundaries into Dirichlet, Neumann, Robin, periodic, symmetry, or open types
  -> Declare interfaces, sheet sources, gauges, and pure-Neumann compatibility
  -> Choose FDM / FVM / FEM / BEM and resolve geometry plus source scales
  -> Assemble operator and inspect symmetry, definiteness, nullspaces, and conservation
  -> Select direct, Krylov, multigrid, transform, or domain-decomposition solver
  -> If charge depends on u: iterate or solve coupled nonlinear residual with damping
  -> Check algebraic residual and integrated source-to-boundary flux balance
  -> Run manufactured solution, analytical limit, and mesh-convergence tests
  -> Evaluate a declared observable with uncertainty and measurement model
  -> Archive equations, boundary map, mesh, tolerances, hashes, and validation evidence

Poisson equation failures can be diagnosed by separating four error layers. A physics error assigns the wrong source, constitutive law, dimension, or boundary behavior. A mathematical error violates compatibility, gauge, interface, or well-posedness. A discretization error misrepresents curvature, flux, geometry, or singularity. An algebraic error stops iterations too early or uses an unsuitable solver. Treating all disagreement by refining the mesh attacks only one layer. The symptom map below keeps remedial action tied to evidence.

One-dimensional Poisson solutions are indispensable unit tests. For $-\kappa u''=f_0$ on $[0,L]$ with constant source and two prescribed endpoint values, the exact solution is a parabola plus a linear term. Its second derivative tests sign and source normalization, its endpoint values test elimination, and its flux difference tests integrated conservation. Layered coefficients give piecewise-linear flux behavior with a derivative jump inversely proportional to $\kappa$ but continuous $\kappa u'$. These cases catch more implementation errors than visually complex demonstrations.

Radial symmetry distinguishes coordinate operators from Cartesian shortcuts. A spherically symmetric field obeys $r^{-2}\partial_r(r^2\partial_ru)$, while cylindrical symmetry uses $r^{-1}\partial_r(r\partial_ru)$. At $r=0$, regularity removes the apparent singularity and implies an even solution with zero radial derivative. Substituting a Cartesian second derivative misses geometric spreading. A uniformly charged sphere yields a quadratic interior potential matched to a $1/r$ exterior, providing a stringent interface and far-field benchmark.

Screened Poisson equations add a reaction scale. The operator $-\nabla\cdot(\kappa\nabla u)+\sigma u=f$ introduces a screening length approximately $\sqrt{\kappa/\sigma}$ when coefficients are uniform and positive. Positive $\sigma$ strengthens coercivity and removes a constant Neumann nullspace, while a negative reaction can make the operator indefinite and permit resonant behavior. Debye–Hückel electrostatics and modified Helmholtz problems share the positive case. Calling every equation with a Laplacian “Poisson” can conceal this altered spectrum and boundary influence.

Anisotropic Poisson operators encode direction-dependent transport. With tensor $\mathbf K$, flux is $-\mathbf K\nabla u$ and need not align with the gradient. Rotating the tensor relative to the mesh introduces cross derivatives and elongated influence. Positive definiteness is required for ellipticity, but a large eigenvalue ratio makes relaxation and interpolation direction sensitive. Line smoothers, semi-coarsening, aligned meshes, or tensor-aware finite elements may be needed. Replacing the tensor with an arithmetic scalar average destroys directional resistance and can shift terminal flux substantially.

Random coefficients turn one solve into an uncertainty ensemble. Spatially uncertain permittivity, conductivity, geometry, doping, or fixed charge makes $u$ and every output random. Monte Carlo is simple but expensive; polynomial chaos, stochastic Galerkin, low-rank, surrogate, and multilevel methods exploit structure under assumptions. Samples must preserve positivity and plausible correlation length. Reporting only the potential from mean inputs generally does not equal the mean potential when coefficients or closure are nonlinear. Solver tolerance should be below sampling and discretization error.

Time-dependent simulations may solve Poisson at every step without making Poisson dynamic. In electroquasistatics, drift–diffusion, particle-in-cell plasma, phase-field, or incompressible flow, evolving charge or constraint data produce a sequence of elliptic solves. Warm starts and reused preconditioners can save work, but accumulated incompatibility or loose solves can violate conservation. The elliptic field responds within the approximation at each step; propagation and displacement-current physics require Maxwell or other dynamic equations when the quasistatic assumption fails.

Pressure Poisson equations enforce incompressibility as a constraint. Projection methods form a source from the divergence of an intermediate velocity, solve for pressure or pressure correction, and subtract its gradient so the updated velocity is divergence free. Pressure boundary conditions derive from momentum and velocity conditions; guessing them can create boundary layers or loss of accuracy. Pure Neumann pressure has an arbitrary constant, which is physically harmless but algebraically must be managed. Variable density produces a variable coefficient and strengthens the analogy with dielectric electrostatics.

Electrostatic force requires an energy- or stress-consistent evaluation. Differentiating field energy with respect to geometry, integrating Maxwell stress on a suitable surface, or integrating charge times field can agree under consistent assumptions. Directly sampling a singular surface field may not. Mesh motion changes both geometry and discrete space, so naive finite differences can include remeshing noise. Self-force from a particle's own deposited charge must be controlled in particle methods. Force convergence can lag potential convergence because it depends on gradients or shape derivatives.

Data interpolation can dominate source error. Doping profiles, tomography, particle clouds, and measured surface charge arrive on grids or locations different from the PDE mesh. Projection should conserve integral source, respect nonnegativity or signed totals as appropriate, and avoid inventing oscillations. Interpolating logarithmic doping as if linear concentration, or averaging a sheet charge into a volume without its thickness, changes the model. The source integral before and after transfer should be recorded as a basic provenance check.

Software benchmarks should test performance and mathematical invariants together. Record assembly time, setup time, solve time, iterations, memory, parallel efficiency, and achieved residual alongside error, balance, and output convergence. A faster solver that stops at a looser effective tolerance is not a fair comparison. Matrix-free methods trade storage for repeated operator work; direct methods trade factor memory for reliable repeated solves. Hardware and library versions matter, but the operator size, polynomial degree, coefficient contrast, and boundary composition are equally important context.

The best formulation is chosen by the hardest feature, not by habit. A regular box with constant coefficient favors transforms or structured finite differences. Complex geometry and heterogeneous media favor finite elements or conservative finite volumes. Infinite homogeneous exteriors favor boundary elements or specialized open boundaries. Strong local refinement favors adaptive unstructured methods. Repeated solves may justify expensive factorization or multigrid setup. The comparison is about conserved quantities, geometry, singularity, coefficient contrast, outputs, and total computational lifecycle rather than which method is universally superior.

The final interpretation should separate potential, field, flux, and source. Potential $u$ is the primary scalar, field is often $-\nabla u$, constitutive flux is often $-\kappa\nabla u$, and source is the divergence of that flux. They have different continuity, units, regularity, and measurement meaning. Potential can remain continuous while field and flux change across an interface; flux can remain continuous while gradient jumps. Plotting all four with consistent sign and units makes interface laws and conservation visible and prevents a smooth potential picture from hiding a wrong field.

Read the Poisson equation through a source-boundary-and-conservation lens rather than a Laplacian-formula-and-potential-plot lens.