Usadel equation of superconductivity

variational formulation
1

Hi all, I’m trying to solve a problem and I am considering if Fenics is a good platform for my needs.

I asked a question on StackExchange (finite element - Solving coupled PDEs with self-consistency condition - Computational Science Stack Exchange) and I’ve been suggested to ask directly in this forum.

I need to code a 2D solver of the Usadel equations of superconductivity. This problem is written for a series of fields f^{s}_n(r) and \vec{f}^{t}_n(r) (which is a set of a finite number, let’s say N, of fields indexed by n) together with a “self-consistent field” \Delta(r). For each n there is a set of coupled PDEs for the f fields

\begin{cases} (D \nabla^2 - 2 \omega_n) f^s_n(r) = - 2 \pi \Delta(r) + 2 i \vec{h}(r)\cdot\vec{f}^t_n(r) \\ (D \nabla^2 - 2 \omega_n) \vec{f}^t_n(r) = 2 i \vec{h}(r) f^s_n(r) \\ \end{cases}

This 2N set of PDEs needs to be solved together with an algebraic non-linear self-consistency equation which looks like

-\Delta(r)\ln(\tau) = 2 \tau \sum_{n=0}^{N} \Big(\frac{\Delta(r)}{2 \omega_n}-f^s_n(r)\Big)

where \vec{h} and \tau are the external parameters.

The usual way to solve the problem is with by finite differences and using a segregated approach: you solve the 2N equations for a given \Delta(r), then you solve the self-consistency equation (e.g. with Newton-Rapson) for \Delta(r) and then you use the result for a new iteration, until convergence.

I would like to solve this problem using finite elements keeping the segregated approach or even as a fully coupled problem. Is it possible or is there some major problem? I get that coding the solver for the PDEs would be fairly easy (they are just PDEs in the end). But is it possible to write easily a solver of the self-consistency equation in Fenics?

Do you know a somewhat similar problem solved in a publicly available code I can look at?