Numerical Methods

This solver allows for direct comparison between several common numerical schemes.

Time Integration Schemes

The time integrator can be selected in SimConfig via the integrator_type variable.

• euler (Forward Euler): A simple, first-order accurate method. It is fast but requires a very small time step to remain stable.
• rk4 (Runge-Kutta 4): A fourth-order accurate method. It is highly accurate and has a much larger stability region than Euler, allowing for significantly larger time steps and faster simulations for a given accuracy.

Spatial Discretization Schemes

The scheme for the parallel gradient can be selected in SimConfig via the scheme_type variable.

Centered Differencing (centered)

• Method: A second-order accurate, non-dissipative scheme.
• Behavior: Excellent at preserving the shape and amplitude of smooth structures. However, it is prone to producing unphysical oscillations (overshoots and undershoots) near sharp gradients, as seen in the comparison below.
• Use Case: Ideal for problems where high accuracy is paramount and the solution is known to be smooth.

Upwind Differencing (upwind)

• Method: A first-order accurate, highly dissipative scheme. The implementation correctly uses a backward difference based on the positive advection velocity.
• Behavior: Very stable and guaranteed not to produce oscillations. This stability comes at the cost of significant numerical diffusion, which smears out sharp features and artificially damps amplitudes.
• Use Case: Useful for problems where stability is the primary concern, or as a component in more advanced schemes.

Blended Differencing (blended)

• Method: A higher-order scheme that adaptively blends the centered and upwind methods. It uses a "smoothness sensor" to detect sharp gradients.
• Behavior: It combines the best of both worlds. It acts like the accurate centered scheme in smooth regions and automatically adds the dissipative upwind scheme near sharp features to suppress oscillations.
• Use Case: The recommended scheme for general-purpose problems, providing a good balance of accuracy and stability.

Interpolation Schemes

The Flux-Coordinate Independent (FCI) method relies heavily on high-quality interpolation. While data is stored on a simple, fixed Cartesian grid (x, y), the physics operators (like the parallel gradient) are calculated in logical flux coordinates (r, θ). The time integrator can be selected in SimConfig via the spline_type variable.

To evaluate the parallel gradient, we need the value of the field at points forward and backward along a magnetic field line. These points almost never coincide exactly with the Cartesian grid nodes. Interpolation is the process of estimating these off-grid values using the surrounding grid data. The choice of interpolation scheme is critical, as it directly impacts the accuracy, stability, and computational cost of the simulation.

Bilinear Interpolation

• Method: The simplest form of multivariate interpolation. It uses a 2x2 stencil (4 grid points) surrounding the query point. The value is first interpolated linearly in one direction (e.g., x) at two locations, and then the final value is found by interpolating linearly between those two results in the other direction (y).
• Accuracy: First-order accurate.
• Properties:
  • Pros: Very fast and simple to implement.
  • Cons: The interpolated function is continuous, but its derivatives are not. This can lead to "kinks" or sharp corners in the gradient field at cell boundaries, introducing significant numerical error and diffusion. It is generally not smooth enough for high-fidelity physics simulations where gradients are important.

Cubic Hermite Splines

• Method: A more advanced form of cubic interpolation that uses a 2x2 stencil. However, unlike other schemes, it requires not only the values at the grid points but also the pre-calculated gradients (∂f/∂x, ∂f/∂y, ∂²f/∂x∂y) at those points.
• Accuracy: Third-order accurate.
• Properties:
  • Pros: Produces a very smooth, C¹-continuous surface (continuous function and first derivatives).
  • Cons: The need to calculate and store gradients at every grid point makes it computationally expensive and memory-intensive. It is often more complex to implement than Catmull-Rom splines for a similar level of smoothness.

Bicubic Interpolation via Catmull-Rom Splines

• Method: This is the scheme used in spline2d.h. It is a powerful form of cubic interpolation that achieves high-quality results using only the function values from a 4x4 stencil (16 grid points). It does not require pre-calculated gradients. It implicitly constructs a smooth curve based on the control points.
• Accuracy: Third-order accurate.
• Properties:
  • Pros:
    • C¹ Continuity: Ensures that the interpolated function and its first derivatives are continuous across cell boundaries. This smoothness is vital for accurately calculating gradients without introducing artificial jumps.
    • Excellent Balance: Offers an ideal trade-off between accuracy, smoothness, and computational cost for physics simulations. It is significantly smoother than Bilinear and less complex to implement than Hermite splines.
  • Cons: Requires a larger stencil than Bilinear, which means boundary condition handling (like the periodic wrapping in this code) is more complex and essential.

Scheme	Stencil Size	Required Data at Nodes	Accuracy	Derivative Continuity
Bilinear	2x2	Values only	First	Discontinuous (C⁰)
Cubic Hermite	2x2	Values and Gradients	Third	Continuous (C¹)
Catmull-Rom (Bicubic)	4x4	Values only	Third	Continuous (C¹)

The choice of Catmull-Rom splines for this solver is deliberate, as it provides the necessary smoothness for calculating physical gradients without the overhead of storing those gradients explicitly.

Boundary Conditions and Performance

• Periodic Boundaries: The 4x4 stencil for Catmull-Rom interpolation requires careful handling at the edges of the domain. The implementation correctly applies periodic boundary conditions by wrapping indices around to the opposite side of the grid when the stencil extends beyond a boundary.
• Optimization: Because interpolation occurs frequently within the main loop, performance is critical. The solver leverages the uniform Cartesian grid to perform fast index lookups using simple arithmetic (index = (x - x_min) / dx), avoiding costly searches.

Comparison of Schemes

The animation below compares the results for the Isotropic Blob problem using the three spatial schemes with an RK4 integrator.

• Centered (Left): Preserves amplitude but shows faint unphysical oscillations around the edges.
• Upwind (Center): No oscillations, but the blob's amplitude rapidly decays due to numerical diffusion.
• Blended (Right): The ideal result. The amplitude is well-preserved, and the oscillations are effectively damped.

---