A question regarding border conditions in coupled PDEs

1

Hello,
as the title says my group partners and I got a question regarding border conditions in Fenics when applied to Coupled PDEs.

At the moment we are working at implementing a simulation of Corona Effects around a high voltage power line.
Since, when we were using the splitting method to implement the connected PDEs, we kept encountering problems with differences between what the handed off data looks like (Matrix of data points) and what the formulas want (vectors of formulas), we thought we could go and try to use the method laid out in the example “ft09_reaction_system.py”.

Sadly that has brought up a question we have been unable to solve so far:
The border Conditions of the two PDEs

bcp_ground = DirichletBC(Q, Constant(0), ground)
bcp_conductor = DirichletBC(Q, Constant(voltage), conductor)
bcp = [bcp_ground, bcp_conductor]

bcd = DirichletBC(Q, Constant(density_constant), conductor)

Till now the first PDE used the Border Condition bcp

phi = TrialFunction(Q)
v = TestFunction(Q)
aP = inner(grad(phi), grad(v))dx
LP = rho_2/epsiv*dx
solve(aP==LP, phi, bcp)

and the second PDE the border condition bcd

rho_2 = TrialFunction(Q)
u = TestFunction(Q)
F = (dot((rho_2**2) / epsi, u) - dot(dot(grad(phi), grad(rho_2)), u))*dx
solve(F==0, rho_2, bcd)

so if we combine those we got one weak form that should look like this

u = TrialFunction(Q)
v = TestFunction(Q)

phi,rho_2 = split(u)
v1,v2 = split(v)

F = inner(grad(phi), grad(v1)) * dx \
+(dot((rho_2**2) / epsi, v2) - dot(dot(grad(phi), grad(rho_2)), v2)) * dx
-rho_2/epsi * v1 * dx

solve(F==0, u, BorderConditionX, Jacobian)

as you have probably spotted there is a slight problem with combining the two PDEs in that they use different Border Conditions which might be defined over the same area but they have definitions on the same border which are different.

So how do we tell fenics to use bcp for one part of the formula or variable and bcd for the other?

On a side note we are also trying to find the right solver for such a highly non linear problem, any tips would be appreciated.

Thank you for your help.

2

If I understand the problem correctly, I believe the solution is to restrict the DirichletBCs bcp and bcd to individual fields within the mixed FunctionSpace, by using Q.sub(0) or Q.sub(1) instead of Q in the definition of the DirichletBCs.

3

Another problem seems to be on the gradient of the phi, normally we should get a vector about direction x and y , but so far our result seems to be a matrix,and I wonder if this has something to do with the rank in the vectorspace or functionspace.

And We have solved phi in the previous step, I wonder if these two lines of code are enough to describe E
E = Function(VectorFunctionSpace)
E.assign(project(-grad(phi), VectorFunctionSpace)

Thanks