Finite volume solution of electrostatics using magnetic vector potential

Question

I would like to solve for the electric potential and magnetic vector potential using the finite volume method (collocated grid). My equations are:

$\nabla\cdot(\sigma\nabla\phi)=0$

$\nabla \cdot \nabla A = \mu_m \sigma \nabla\phi$

The domain is a wedge, with the following boundary conditions:

• A fixed gradient for the electric potential $\phi$ on the top, zero gradient at the bottom and a zero value at the side boundary
• Zero gradient for the magnetic vector potential $A$ on all boundaries

I didn't expect this would be a tough problem, but the linear solver (GAMG) only converges with under-relaxation, and the result is dependent on the value of the under-relaxation factor. Can anyone give me a hint where things are going wrong?

Additional information: The software package I use is OpenFOAM. The gradients are discretized using the Gauss linear scheme, and the laplacian with "Gauss linear corrected" (explicit non-orthogonal correction).

Answer 1

The solution required two changes:

• Change the boundary conditions for $A$ to fixedValue 0, which follows from $B\cdot n = 0$ at the boundaries.
• Disable under-relaxation. Even though the system converged without under-relaxation, I was not able to reach convergence no matter which value I assigned to the under-relaxation coefficient. This seems to be a problem with the implementation of matrix relaxation in OpenFOAM.

Answer 2

OpenFOAM tends to follow [Ferziger & Peric, Computational Methods for Fluid Dynamics, Springer, 2002] so my guess is that they solve

\begin{align} & \text{outer iter: do $k=0,K-1$} \\ & \qquad \text{inner iter: solve} \quad \left( A_k + \frac{1-\omega}{\omega} \mathrm{diag}(A_k) \right) x_{k+1} = b + \frac{1-\omega}{\omega} \mathrm{diag}(A_k) x_k \\ & \text{end outer iter} \end{align}

Upon outer-loop convergence $x_{k+1}=x_k$ so the correct system is being solved. But in your linear case, you don't have an outer iter (k=0 with initial guess $x_0=0$) so you get the solution of

\begin{equation} \left( A + \frac{1-\omega}{\omega} \mathrm{diag}(A) \right) x = b \end{equation}

which is only correct for relaxation parameter $\omega=1$.

• To confirm my guess, try to construct $A$ and substract $\frac{1-\omega}{\omega} \mathrm{diag}(A)$ before passing it to the linear solver.

• To actually use their relaxation scheme, try to add the outer iterations (take for example K=10 and solve the linear system with a relative tolerance of 1e-2 instead of taking K=1 and solving with a relative tolerance of 1e-12)