This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

Solve domains sequentially: how to estimate the BC at the interface

0 votes

Let's assume we have two domains A and B. Their interface is noted I.
One Dirichlet BC is applied at the left part of A, one other Dirichlet BC is applied at the right part of B (stationary diffusion problem)
If you consider the whole domain A+B, you can solve right now the problem, there is no issue (you need only to mark your two sub domains and then allocate the diffusion coefficient accordingly)

Now, I want to solve the diffusion in both domains sequentially: first solve domain A, then solve domain B (with 2 distinct variational formulations: one for each domain) and iterate until convergence. Do not ask why, it's just a test case i try to solve for a much more complex problem.
Here, the issue is you do not have the BC at the interface I.
Instead of, we know the concentration flux should be equal at the interface I.

My approach is:
- set a Dirichlet BC at the interface (initial guess), identical for both domain
loop until residual is 0:
- Solve domain A with current Dirichlet BC at the interface
- Solve domain B with current Dirichlet BC at the interface
- Calculate the flux discontinuity at the interface (that's the residual i want to minimize)
- Estimate new values for the Dirichlet BC at the interface
- current Dirichlet BC = new values for the Dirichlet BC at the interface

Then, the function to minimize (evaluated at the interface only) is
F = - Da * grad(Ca).nia + Db * grad(Cb).nib
with Dx the diffusion coefficient of domain x
Cx the concentration of domain x
nix the unit normal vector at the interface, for the domain x

So, if the interface is flat, then the problem is actually 1D, and you can use a newton method even without knowing F (you calculate the flux discontinuity just before and after the initial guess to get a numerical derivate and that's all): dx = x1-x0 = -f(x0)/f'(x0)

But, i do not have a flat interface, then instead of i should calculate the Jacobian of F and solve Jacobian(F(x0)).dx = -F(x0), such as:

A,b = assemble_system(J, -F, bc0) 
solve(A,dx.vector(), b)

I have difficulties to define correctly J, F and bco
Some help please?
Regards,

asked Nov 10, 2016 by fussegli FEniCS Novice (700 points)

Hi, I have solved my problem by calculating term by term the Jacobian matrix...
It works, but it's definitely not the smartest way to do it I am sure.
If someone has an idea how to do it.
Thx.

...