How to solve coupled PDEs for electrostatics?



I am looking into solving the 2 following coupled PDEs:
(i) Div[Grad(V)] = -q/cste
(ii) Div[q.Grad(V)] = 0
With q and V being functions: q(x;y) & V(x;y)
With boundary conditions on V all around the 2D domain

  • (i) is a basic Poisson problem and (ii) is a variable coefficient Poisson problem
    Independant resolution of (i) and (ii) works fine with FEniCS, thanks to relevant documentation.
    For both cases: Assuming q is known as a starting point, I can get V as a result on my mesh/domain
  • I tried to couple (i) & (ii) as follows
    STEP1: Assume a first q, set BCs on V and solve (ii) => V1 as an output
    STEP2: Using V1, calculate q with (i) => q2 as an output
    STEP3: Using q2, set (same) BCs on V and solve (ii) => V3 as an output
    Unfortunately, result doesn’t converge
  • I also tried to work with some coupling FEniCS tutorials, but they are more complicated than my need and I couldn’t adapt them

Does anyone know a simple way to couple (i) and (ii) in the same variational statement, and get q and v efficiently ?

In advance Thank you for your help !