Skip to content

Electric Solver Module

Overview

The electric solver module provides numerical solutions for electrostatic field problems within the OpenPhase framework. It is designed to compute electric potentials and derived quantities (e.g. electric field, current density) on a structured grid, consistent with the global simulation domain.

The module follows an abstract interface (ElectricSolver) and currently includes a spectral implementation (ElectricSolverSpectral) based on FFTW.

This document explains the architecture, numerical background, and extension points. It is not a full API reference.


Mathematical Background

The solver addresses electrostatic problems of the form:

(1)(σ(x)ϕ(x))=0

where:

  • ϕ — electric potential
  • σ — electric conductivity (possibly spatially varying)
  • E=ϕ — electric field

For homogeneous conductivity and periodic boundary conditions, this reduces to a Poisson-type equation that can be efficiently solved in Fourier space.

The spectral solver leverages the identity:

(2)2ϕF|k|2ϕ^

which allows solving the equation algebraically in frequency space.


Key Classes and Concepts

  • ElectricSolver (abstract base) — defines the common solver interface: virtual void Solve(ElectricProperties& EP, BoundaryConditions& BC) = 0;. To add a new solver, derive from this class and implement Solve().
  • ElectricSolverSpectral — the FFT-based concrete implementation. Uses FFTW3 for forward / backward transforms on structured Cartesian grids; optimised for periodic boundaries and large 3D domains. Constructed as ElectricSolverSpectral(Settings& locSettings, const std::string InputFileName), which allocates the real- and Fourier-space arrays, creates the FFTW plans, and extracts the grid dimensions from Settings.
  • ElectricalPotential — owns the electrode list and per-component elementary charges. See Electrics/ElectricalPotential.md.
  • ElectricProperties — per-component charge and polarization density storage. See Electrics/ElectricProperties.md.

Internal Data Layout

The solver uses flattened 3D indexing:

cpp
size_t rlidx(long int i, long int j, long int k);

Mapping:

[ XYZ = k + Nz (j + Ny i) ]

This layout ensures:

  • Cache-friendly memory access
  • Compatibility with FFTW real-to-complex transforms
  • Direct mapping between physical and spectral domains

Memory Management

Internal arrays:

  • rhs — real-space right-hand side
  • freq — complex Fourier-space buffer
  • ForwardPlan — FFTW forward transform
  • BackwardPlan — FFTW backward transform

Memory is allocated during Initialize() and freed in the destructor.

When extending the solver, ensure:

  • FFTW plans are destroyed properly
  • Arrays are deallocated before reallocation
  • No grid-size mismatch occurs after remeshing

Solver Workflow

The Solve() method follows this high-level procedure:

  1. Assemble real-space RHS (based on conductivity and boundary conditions)
  2. Execute forward FFT
  3. Solve algebraic system in Fourier space
  4. Execute inverse FFT
  5. Normalize result
  6. Update electric properties (potential, field)

Pseudo-algorithm:

text
Build RHS in real space
FFT → frequency space
For each k ≠ 0:
    φ̂(k) = RHŜ(k) / (-|k|²)
Handle k = 0 mode (reference potential)
Inverse FFT
Normalize
Compute electric field

Boundary Conditions

The spectral solver naturally supports:

  • Fully periodic boundary conditions

Non-periodic boundary conditions require either:

  • Modification of the RHS construction
  • Green's function approaches
  • Alternative solvers (e.g., finite difference or FEM)

When implementing new boundary types, consider deriving a separate solver class.


Grid Parameters

The solver uses GridParameters:

  • Nx, Ny, Nz
  • Physical dimensions
  • Grid spacing

Consistency with the global simulation grid is mandatory.

Changing grid size requires:

  • Reallocation of FFT buffers
  • Recreation of FFT plans

Extension Guidelines

Implementing a New Electric Solver

  1. Derive from ElectricSolver

  2. Implement:

    cpp
    void Solve(ElectricProperties& EP,
               BoundaryConditions& BC) override;
  3. Register and initialize in the same manner as other OPObjects

  4. Ensure compatibility with:

    • ElectricProperties
    • BoundaryConditions
    • Settings

When to Use the Spectral Solver

Recommended for:

  • Periodic domains
  • Large 3D grids
  • Homogeneous or smoothly varying conductivity
  • HPC environments

Not ideal for:

  • Strongly localized boundary conditions
  • Complex non-periodic geometries
  • Highly anisotropic conductivity tensors (unless extended)

Performance Considerations

The spectral solver scales as:

(3)O(NlogN)

where N=NxNyNz.

Performance tips:

  • Use power-of-two grid sizes where possible
  • Reuse FFTW plans
  • Avoid frequent reinitialization
  • Compile FFTW with optimization flags
  • Enable FFTW wisdom for repeated runs

Numerical Considerations

  • The zero-frequency mode must be handled carefully.
  • Potential is defined up to an additive constant.
  • For charge-neutral systems, enforce mean-zero potential.
  • Aliasing effects can occur if the RHS contains sharp discontinuities.

Usage

Input

ElectricSolverSpectral has no dedicated .opi block. The solver reads its inputs from the other electrics classes:

  • Per-component elementary-charge valence from ElectricalPotential — the $Nelectrodes, $CntElectr<i>, $Z<q> block.
  • Charge density from ElectricProperties, computed via CalculateChargeDensity(Composition&).
  • Face-level boundary handling from BoundaryConditions — periodic is the natural fit for the spectral solver.

Output

  • VTK fields for the electric potential and derived electric field through WriteVTK.
  • Field diagnostics through console messages (ConsoleOutput).

Example

cpp
#include "ElectricalPotential.h"
#include "Electrics/ElectricProperties.h"
#include "Electrics/ElectricSolverSpectral.h"

ElectricalPotential    EPot(OPSettings, InputFile);
ElectricProperties     ElecP(OPSettings, InputFile);
ElectricSolverSpectral ElecS(OPSettings, InputFile);

for(RTC.tStep = RTC.tStart; RTC.tStep <= RTC.nSteps; RTC.IncrementTimeStep())
{
    ElecP.CalculateChargeDensity(Cx);
    ElecP.SetBoundaryConditions(BC);
    ElecS.Solve(ElecP, BC);

    if (RTC.WriteVTK())
    {
        ElecS.WriteVTK(OPSettings, RTC.tStep);
    }
}

Dependencies

Summary

The electric solver module provides:

  • A clean abstract interface (ElectricSolver)
  • A high-performance FFT-based implementation (ElectricSolverSpectral)
  • Modular design compatible with the OpenPhase architecture
  • Efficient large-scale electrostatic field computation

It is designed for extensibility and high-performance scientific computing environments.

Released under the GNU GPLv3 License.