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

4 months ago by
Hello there,
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$κ1T1c1u·T1=0                      in  $\Omega_1$Ω1   
$\nabla\kappa_1\nabla T_1-c_1u\cdot\nabla T_1=h\left(T_2-T_1\right)$κ1T1c1u·T1=h(T2T1)     in  $\Omega_2$Ω2
$\nabla\kappa_2\nabla T_2=Q+h\left(T_2-T1\right)$κ2T2=Q+h(T2T1)                  in $\Omega_2$Ω2
where  $\kappa_1$κ1 ,  $\kappa_2$κ2 ,  $c_1$c1 , 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)

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

Community: FEniCS Project
What I tried was that, in  $\Omega_1$Ω1  as I added a new PDE
$\nabla\kappa_f\nabla T_2=0$κƒ T2=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}$κƒ =1e10
Now I can solve the PDE, but how big difference is compared with the original problem?
written 4 months ago by SICHENG SUN  
This recent question might have a couple of links (disclaimer: one of them is mine). The development version dolfinx will also probably address this issue in the next future.
written 3 months ago by Francesco Ballarin  
I am sorry I did not see your reply until now. Thank you so much. I will take a look into it.
written 3 months ago by SICHENG SUN  

1 Answer

3 months ago by
Bio Pi  
Hi there, I have a similar problem. Just to be clear,  $\Omega_1\subset\Omega_2$Ω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,


  $D_1\nabla^2u_1+f\left(u_2\right)=0\text{ in }\Omega_1$D12u1+ƒ (u2)=0 in Ω1  
$u_1=u^{boundary}_1\text{ in }\partial\Omega_1$u1=uboundary1 in Ω1
$u_1=0\text{ otherwise }$u1=0 otherwise
$D_2\nabla^2u_2+g\left(u_1\right)=0\text{ in }\Omega_2$D22u2+g(u1)=0 in Ω2
$u_2=u^{boundary}_2\text{ in }\partial\Omega_2$u2=uboundary2 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$uboundary1=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$D1=0 and ƒ 1(u2)=0 outside of Ω1 . This solution is not general and nor elegant nor efficient AND  more importantly I dont know if is correct.

Hello, If you are still working on this problem, I think you may want to check ".ident_zeros()", there is a relative thread:
However, I myself have not yet figure out how to use this command.
written 3 months ago by SICHENG SUN  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »