I don’t know where to begin, I have been working on a single problem since several years now, and it works well in a certain regime (which I want to extend now), but I am facing numerous difficulties. I am at the point of wondering whether I should change my approach and use RT elements and introducing a new unknown.
In short: I solve for T and V in the equations \nabla \cdot \vec J_U=0 and \nabla \cdot \vec J=0 where \vec J=−\sigma \nabla V−\sigma S \nabla T and \vec J_U = -\kappa \nabla T +TS\vec J +V\vec J. My weak form is:
# Weak form.
J_vector = -sigma * grad(volt) - sigma * Seebeck_tensor * grad(temp)
Fourier_term = dot(-κ * grad(temp), grad(u)) * dx
F_T = Fourier_term
rest_terms = dot(temp * Seebeck_tensor* J_vector, grad(u)) * dx + dot(volt * J_vector, grad(u)) * dx
F_T += rest_terms
F_V = -dot(grad(v), sigma * grad(volt))*dx -dot(grad(v), sigma * Seebeck_tensor * grad(temp))*dx
weak_form = F_T + F_V
which corresponds to \nabla \cdot (\kappa \nabla T ) +\nabla \cdot (T S \vec J)+\nabla \cdot(V\vec J)=0 and \nabla \cdot \vec J=0, where I set Dirichlet boundary conditions for T on 2 boundaries (out of 4) and potentially a Dirichlet b.c. at a single point for V (this does help a bit for Newton method to converge).
Here is the thing, S is a tensor, a 2x2 diagonal matrix with entries S_{xx} and S_{xy}. And it turns out that in the regime where this anisotropy is “small” enough not to perturb T much, then my FEniCSx code runs fine and is able to retrieve some analytical formulae I had derived. The problem arises when the anisotropy in S becomes big enough so that T starts to become impacted. For example, in my code I compute:
assemble_scalar(form(div(J_U) * dx))
which should be 0 since this corresponds to the heat equation I am solving. Well, this number starts to increase (it is about 1e-7 or even 1e-8 in the regime of a small anisotropy, where my code runs fine), it becomes something like -0.03 when anisotropy is big and then Newton method diverges. I did implement a trick to make Newton method converge, by solving many intermediate problems with increasing anisotropy in S, but this doesn’t change noticeably the result, i.e. it helps at converging, but it converges to the same bogus values. (I know I am not computing the L2 norm, maybe I should do it?).
As regards to the volume integral of \nabla \cdot \vec J, it also starts to depart from 0 but it is closer to 0 than its \vec J_U counterpart is.
I do not know whether my solution for T and V is accurate, or if the inaccuracy is due to computing gradients.
What I know is that I do need to compute the volume integral of \vec J^2 and spatial derivatives of \vec J, again, integrated over the whole volume. So, I would need an accurate \vec J, I suppose.
I am willing to introduce new variables such as J, J_U, J_Q (some heat flux) and use RT elements using a ME scheme if it makes things better. I tested rapidly and I got even worse results by implementing J with RT elements, where the div J and div J_U values where a million times greater than in my case, for the regime where my code works fine.
Any pointer or idea is welcome. My full code is gigantic and a MWE cannot be implemented, as these values are only achieved in my specific problem. I don’t mind to post my full code here, it’s just big, and the mesh comes from a (g)msh file.
