A nonlinear Harmonic Navier-Stokes (HNS) framework is introduced for simulating instabilities in laminar spanwise-invariant shear layers, featuring sharp and smooth wall surface protuberances. While such cases play a critical role in the process of laminar-to-turbulent transition, classical stability theory analyses such as parabolized or local stability methods fail to provide (accurate) results, due to their underlying assumptions. The generalized incompressible Navier-Stokes (NS) equations are expanded in perturbed form, using a spanwise and temporal Fourier ansatz for flow perturbations. The resulting equations are discretized using spectral collocation in the wall-normal direction and finite-difference methods in the streamwise direction. The equations are then solved using a direct sparse-matrix solver. The nonlinear mode interaction terms are converged iteratively. The solution implementation makes use of a generalized domain transformation to account for geometrical smooth surface features, such as humps. No-slip conditions can be embedded in the interior domain to account for the presence of sharp surface features such as forward- or backward-facing steps. Common difficulties with Navier-Stokes solvers, such as the treatment of the outflow boundary and convergence of nonlinear terms, are considered in detail. The performance of the developed solver is evaluated against several cases of representative boundary layer instability growth, including linear and nonlinear growth of Tollmien-Schlichting waves in a Blasius boundary layer and stationary crossflow instabilities in a swept flat-plate boundary layer. The latter problem is also treated in the presence of a geometrical smooth hump and a sharp forward-facing step at the wall. HNS simulation results, such as perturbation amplitudes, growth rates, and shape functions, are compared to benchmark flow stability analysis methods such as Parabolized Stability Equations (PSE), Adaptive Harmonic Linearized Navier-Stokes (AHLNS), or Direct Numerical Simulations (DNS). Good agreement is observed in all cases. The HNS solver is subjected to a grid convergence study and a simple performance benchmark, namely memory usage and computational cost. The computational cost is found to be considerably lower than high-fidelity DNS at comparable grid resolutions. Program summaryProgram Title: DeHNSSoCPC Library link to program files:https://doi.org/10.17632/9bnms99kk2.1Developer's repository link:https://github.com/SvenWesterbeek/DeHNSSoLicensing provisions: GPLv3Programming language: MatlabSupplementary material: The supplementary material contains the code as well as a user manual.Nature of problem: Fluid flows are subject to laminar-to-turbulent transition following the growth of instabilities. To avoid computationally demanding Direct Numerical Simulations (DNS), perturbation theory is often applied to their analysis. However, classical stability methods based on the Orr-Sommerfeld equation or the Parabolized Stability Equations neglect the influence of streamwise gradients in varying degrees. The validity of these assumptions is difficult to estimate a priori.Solution method: The Delft Harmonic Navier-Stokes Solver (DeHNSSo) solves the harmonic Navier-Stokes equations nonlinearly on domains featuring sharp and smooth spanwise invariant surface features using a generalized grid approach in combination with an embedded boundary method. This allows the user to include the effects of streamwise gradients on flow at a fraction of the cost of DNS.Additional comments including restrictions and unusual features: In DeHNSSo, the equations are solved using direct matrix solvers. As such, memory is treated as a scarce resource. The problem is formulated in a mode-independent manner such that left-hand side matrices need only be computed and stored once despite incorporating nonlinear terms. Additionally, DeHNSSo offers the user the possibility to prescribe inhomogeneous boundary conditions to introduce and solve the receptivity problem or manipulate instabilities. Due to the double Fourier expansion, the solver is restricted to spanwise and temporally periodic problems.