### How to solve coupled PDEs with one euqation only defined on a submesh

I am facing a problem which requires to solve a coupled PDE, but one of the equation is only defined on the submesh.

The problem is stated below:

$\nabla\kappa_1\nabla T_1-c_1u\cdot\nabla T_1=0$∇

`κ`

_{1}∇

`T`

_{1}−

`c`

_{1}

`u`·∇

`T`

_{1}=0 in $\Omega_1$Ω

_{1}

$\nabla\kappa_1\nabla T_1-c_1u\cdot\nabla T_1=h\left(T_2-T_1\right)$∇

`κ`

_{1}∇

`T`

_{1}−

`c`

_{1}

`u`·∇

`T`

_{1}=

`h`(

`T`

_{2}−

`T`

_{1}) in $\Omega_2$Ω

_{2}

$\nabla\kappa_2\nabla T_2=Q+h\left(T_2-T1\right)$∇

`κ`

_{2}∇

`T`

_{2}=

`Q`+

`h`(

`T`

_{2}−

`T`1) in $\Omega_2$Ω

_{2}

where $\kappa_1$

`κ`

_{1}, $\kappa_2$

`κ`

_{2}, $c_1$

`c`

_{1}, h are just some coefficient. u is the velocity I solved from a Navier Stokes equations. T_1 and T_2 here are my variable.

Thus in my fenics code, I have:

```
Q_ele = FiniteElement('CG', mesh.ufl_cell(), 1)
ME_ele = Q_ele*Q_ele
ME = FunctionSpace(mesh, ME_ele)
s1, s2 = TestFunctions(ME)
T = Function(ME)
T1, T2 = split(T)
F = k1*inner(grad(T1), grad(s1))*dx(1) + c1*inner(u, grad(T1))*s1*dx(1)\
+ k1*inner(grad(T1), grad(s1))*dx(2) + c1*inner(u, grad(T1))*s1*dx(2)\
- h*(T2-T1)*s1*dx(2)\ # weak form of 2nd pde
+ k2*inner(grad(T2), grad(s2))*dx(2)\
+ h*(T2-T1)*s2*dx(2) + Q*s2*dx(2) # weak form of my 3rd pde
J = derivative(F, T)
problem = NonlinearVariationalProblem(F, T, bc, J)
problem_solver = NonlinearVariationalSolver(problem)
problem_solver.solve()
```

But the newton iteration blows up at even the 1st step with the residual goes to -nan.

I guess the reason is that my second equation with respect to T_2 has no contribution in $\Omega_1$Ω_{1} , but T_2 is defined in the whole domain. Any suggestions about how to solve it?

Thank you

### 1 Answer

_{1}⊂Ω

_{2}, right? Also, if you say that $\kappa$

`κ`'s are only coeficients, what the gradients means? In your case are zero and your equations simplify. I would start with a toy problem such as,

and

$D_1\nabla^2u_1+f\left(u_2\right)=0\text{ in }\Omega_1$`D`_{1}∇^{2}`u`_{1}+`ƒ` (`u`_{2})=0 in Ω_{1}

$u_1=u^{boundary}_1\text{ in }\partial\Omega_1$`u`_{1}=`u`^{boundary}_{1} in ∂Ω_{1}

$u_1=0\text{ otherwise }$`u`_{1}=0 otherwise

and

$D_2\nabla^2u_2+g\left(u_1\right)=0\text{ in }\Omega_2$`D`_{2}∇^{2}`u`_{2}+`g`(`u`_{1})=0 in Ω_{2}

$u_2=u^{boundary}_2\text{ in }\partial\Omega_2$`u`_{2}=`u`^{boundary}_{2} in ∂Ω_{2}

you can try to solve this problem in 1D and check if your code is right.

In my case $u^{boundary}_1=0$`u`^{boundary}_{1}=0, so I partialy "workaround" this issue by introducing $D_1=0\text{ and }f_1\left(u_2\right)=0\text{ outside of }\Omega_1$`D`_{1}=0 and `ƒ` _{1}(`u`_{2})=0 outside of Ω_{1} . This solution is not general and nor elegant nor efficient AND more importantly I dont know if is correct.

https://answers.launchpad.net/fenics/+question/188087

However, I myself have not yet figure out how to use this command.

_{1}as I added a new PDE$\nabla\kappa_f\nabla T_2=0$∇

κ_{ƒ }∇T_{2}=0 , but this will change the T_2 in $\Omega_2$Ω_{2}as some heat may flow from Ω2 to Ω1.Thus in order to prevent the heat flux from Ω

_{2}to Ω_{1 }, I choose a very small $\kappa_f=1e^{-10}$κ_{ƒ }=1e^{−10}Now I can solve the PDE, but how big difference is compared with the original problem?