Deriving the weak form of my problem, is it correct?

[raboynics]

I am solving 2 coupled PDEs in a single material with anisotropic physical property S =\begin{bmatrix} S_{xx} & 0 \\ 0 & S_{yy}\end{bmatrix}, all other physical properties are scalar constants. The region being a quarter of an annulus, with therefore 4 boundaries.

It is question of the heat equation \nabla \cdot \vec J_U=0 and the conservation of electrical charge in steady state: \nabla \cdot \vec J =0.

Where \vec J_U =-\kappa \nabla T + ST\vec J + V\vec J and \vec J =-\sigma \nabla V - \sigma S \nabla T. I am therefore solving for V and T with finite elements with a mixed element scheme. My goal is to build \vec J and \vec J_U once V and T are solved for, so that I can compute these fluxes along the boundaries.

The boundary conditions are a fixed temperature on 2 boundaries, i.e. Dirichlet boundary conditions for T on 2 sides. The material being thermally isolated, I also want natural Neumann boundary conditions for T on the remaining 2 sides, and I want no entering electric current on any side, which means \vec J has to be parallel to all boundaries. This translates as a complicated value of V on the boundaries, but which aren’t enforced as neither Neumann nor Dirichlet b.c.s are far as I know. FEniCSx seems to be able to figure it out pretty well with the following weak form:

J_vector = -sigma * grad(volt) - sigma * S_tensor * grad(temp)
F_T = dot(-κ * grad(temp), grad(u)) * dx + dot(sigma * temp * S_tensor* J_vector, grad(u)) * dx + dot(volt * J_vector, grad(u)) * dx
F_V = -dot(grad(v), sigma * grad(volt))*dx -dot(grad(v), sigma * S_tensor * grad(temp))*dx
weak_form = F_T + F_V

which I reached by doing \nabla \cdot \vec J_U = \nabla \cdot (-\kappa \nabla T) +\nabla \cdot (TS\vec J) + \nabla \cdot (V\vec J)=0, integrated by parts the 2nd and 3rd terms, assumed that the “test function” u (which isn’t really a test function in FEniCSx since I use a mixed element space) vanishes on the boundaries. This means I assumed \int _{\partial \Omega}uST\vec J \cdot \hat n ds = \int _{\partial \Omega}uV\vec J \cdot \hat n ds=0. I am not sure whether I can assume this… Hence me posting here.
And F_V is found from the 2nd PDE’s with a simple integration by parts.

I obtain a result that makes a lot of sense to me, it agrees pretty well to an approximate solution to this problem (but without the coupling between V and T). I believe the implementation is correct, but I am not 100% sure because of the fluxes I may have ignored by assuming the test functions to vanish on the boundaries.

Also, I do not quite understand what is going on when I do not specify any Neumann boundary conditions. “Natural” ones are assumes, yes, okay, but what does it mean here? Which flux(es) are zero on the boundaries if I do not provide any Neumann b.c.s ?

For information, the obtained \vec J is displayed on the picture below. It agrees almost perfectly with the expected one.

Here is V and some isolines. It seems like The isolines are perpendicular to the boundaries where no thermal boundary conditions have been applied.

And for completeness, here is T, so it’s pretty clear which boundaries have Dirichlet boundary conditions applied (the inner and outer radii).

Edit: I think I start to understand what happens regarding which fluxes are zero if I don’t specify any boundary condition (neither Dirichlet nor Neumann). Here, since the gradient is perpendicular to isolines, it is evident that on the 2 boundaries without any boundary conditions, both \nabla T and \nabla V are parallel to the boundaries, meaning they have no normal component there.

The case of V is a bit more special than that of T because V has absolutely no boundary condition specified at all. The gradient of V is oblique to the curved boundaries where T is kept fixed, whereas \nabla T is perpendicular to these boundaries, which looks good to me.

[Stein]

@raboynics

F_T = dot(-κ * grad(temp), grad(u)) * dx + dot(sigma * temp * S_tensor* J_....

The sigma seems incorrect in that second integral (is not in your strong form)

@raboynics

integrated by parts the 2nd and 3rd terms

You did integration by parts on all three terms on the right-hand side.

@raboynics

assumed that the “test function” u (

Here you have some major misconceptions about weak formulations and FEM. Your testfunction vanishes only where you specify a Dirichlet boundary condition on the associated trial space. If you did not do this, then you cannot assume it vanishes (it simply does not). Which leads to…

@raboynics

This means I assumed \int _{\partial \Omega}uST\vec J \cdot \hat n ds = \int _{\partial \Omega}uV\vec J \cdot \hat n ds=0. I

Per my post in the other thread, this is not an assumption. This is your imposition of the natural boundary condition. So not u = 0, but the normal component of the flux is zero \vec J_U \cdot \hat n= 0 and \vec J \cdot \hat n= 0. That is a fundamental concept in weak formulations of PDEs.

You seem to be making it more complicated in your head than it is. Say you are solving:
\nabla \cdot J = 0
Multiplication with test, integration over domain, and integration by parts gives:
\int_\Omega J\cdot\nabla v = \int_{\partial\Omega} J\cdot n v
So if you are actually solving:
\int_\Omega J\cdot\nabla v = 0
then this is implying:
\int_{\partial\Omega} J\cdot n v = 0
If you did not specifiy Dirichlet conditions (which would mean v=0 there), then this means that you are imposing J\cdot n = 0,

[raboynics]

@Stein

The sigma seems incorrect in that second integral (is not in your strong form)

Oops, thanks a lot for pointing this out. I didn’t notice it because sigma is worth 1 by default in my code, but I am currently working on a region-dependent sigma, so this will matter a lot.

@Stein

You did integration by parts on all three terms on the right-hand side.

You’re right. In my last derivations on my draft I only focused on the last 2 terms, and I had forgotten I actually did the same for the first term also, in the past.

@Stein

Here you have some major misconceptions about weak formulations and FEM. Your testfunction vanishes only where you specify a Dirichlet boundary condition on the associated trial space. If you did not do this, then you cannot assume it vanishes (it simply does not). Which leads to…

Alright.

@Stein

Per my post in the other thread, this is not an assumption. This is your imposition of the natural boundary condition. So not u = 0, but the normal component of the flux is zero \vec J_U \cdot \hat n= 0 and \vec J \cdot \hat n= 0. That is a fundamental concept in weak formulations of PDEs.

You seem to be making it more complicated in your head than it is. Say you are solving:
\nabla \cdot J = 0
Multiplication with test, integration over domain, and integration by parts gives:
\int_\Omega J\cdot\nabla v = \int_{\partial\Omega} J\cdot n v
So if you are actually solving:
\int_\Omega J\cdot\nabla v = 0
then this is implying:
\int_{\partial\Omega} J\cdot n v = 0
If you did not specifiy Dirichlet conditions (which would mean v=0 there), then this means that you are imposing J\cdot n = 0,

Ok… So here, I must think quite a bit. I do want \vec J to be parallel to all boundaries, therefore I am 100% satisfied with imposing \vec J \cdot \hat n=0 on all boundaries. This corresponds to what happens physically.

However things are more complicated regarding \vec J_U. I think something strange happens physically (but I am not 100% sure yet), if I apply Dirichlet boundary conditions on T at the \theta=0 and \theta=\pi/2 boundaries. In that case, \vec J is still parallel to all boundaries, \vec J_U has a normal component on both the \theta=0 and \theta=\pi/2 boundaries, but I also believe (not yet 100% sure) it should possess a non zero normal component on the r=r_i and r=r_o boundaries where no Dirichlet boundary conditions are applied, for neither T nor V. Therefore I believe \vec J_U \cdot \hat n \neq 0 on these 2 boundaries even though no Dirichlet boundary conditions have been specified. And I do believe FEniCSx points that this indeed holds.

I am not sure if I am being clear…

But I may be missing the surface or boundaries integrals \int_{\partial \Omega}TS\vec J \cdot \hat n ds and \int_{\partial \Omega}V\vec J \cdot \hat n ds. Maybe I should add them in the weak form, indeed.

But then… I do not understand how FEniCSx is able to figure out that the STJ flux term does not vanish on those 2 boundaries where I didn’t prescribe any Dirichlet b.c.s. So yeah, I am still not understanding this part well.

EDIT: I think I understand. I just got “lucky” in that FEniCSx applied the correct conditions. Let me explain.

For some reason FEniCSx correctly applied the \vec J \cdot \hat n =0 condition on all boundaries, just as I wanted and just as I specified in my problem. Good.

However, by neglecting the 2 boundaries integrals \int_{\partial \Omega}TS\vec J \cdot \hat n ds and \int_{\partial \Omega}V\vec J \cdot \hat n ds, I am telling FEniCSx I want TS\vec J \cdot \hat n =0 and V\vec J \hat n =0 on the 2 boundaries with no Dirichlet b.c.s. I do not particularly want to impose this, so this is bad! But… it turns out that this condition is not compatible with the first condition, the one on the normal component of \vec J. Indeed, if \vec J is parallel to the boundaries as I want and as I specified, then S\vec J will NOT be parallel to the boundaries, and so the condition TS\vec J \cdot \hat n =0 is not possible to satisfy.

It turns out I probably got lucky and FEniCSx, for some reason, chose to satisfy the \vec J \cdot \hat n =0 condition over the TS\vec J \cdot \hat n =0 condition (which I would not want to apply!).

What do you think @Stein? I hope my understanding improved

[Stein]

@raboynics

But I may be missing the surface or boundaries integrals \int_{\partial \Omega}TS\vec J \cdot \hat n ds and \int_{\partial \Omega}V\vec J \cdot \hat n ds. Maybe I should add them in the weak form, indeed.

Oh, no, you wouldn’t be ‘allowed’ to put those in the weak form, that would also be a variational crime (just like you can’t operate the Laplacian on an H1 function, you are equally not permitted to evaluate a gradient of such a function on the boundary). That’s exactly why you have to replace those via Natural or Essential BCs.

@raboynics

the \vec J \cdot \hat n =0 condition on all boundaries, just as I wanted and just as I specified in my problem. Good.

However, by neglecting the 2 boundaries integrals \int_{\partial \Omega}TS\vec J \cdot \hat n ds and \int_{\partial \Omega}V\vec J \cdot \hat n ds, I am telling FEniCSx I want TS\vec J \cdot \hat n =0 and V\vec J \hat n =0 on the 2 boundaries with no Dirichlet b.c.s.

You’re confusing me with all the different fluxes again I wrote J\cdot n = 0 and J_U \cdot n =0. That’s two conditions. You’re writing J\cdot n =0 and TSJ\cdot n= 0 and V J \cdot n=0 . I don’t understand where those come from. (And they are all zero if J\cdot n = 0 btw).

Just to clarify; the entire boundary integral you obtain from integration by parts is zero. You can’t say anything about each component of the sum individually.

[raboynics]

@Stein

Oh, no, you wouldn’t be ‘allowed’ to put those in the weak form, that would also be a variational crime (just like you can’t operate the Laplacian on an H1 function, you are equally not permitted to evaluate a gradient of such a function on the boundary). That’s exactly why you have to replace those via Natural or Essential BCs.

Ok… so, I am sorry but I do not understand how to correctly derive the weak form I seek. At first glance though, it looks like FEniCSx nails it perfectly well though.

@Stein

You’re confusing me with all the different fluxes again I wrote J\cdot n = 0 and J_U \cdot n =0. That’s two conditions. You’re writing J\cdot n =0 and TSJ\cdot n= 0 and V J \cdot n=0 . I don’t understand where those come from. (And they are all zero if J\cdot n = 0 btw).

Just to clarify; the entire boundary integral you obtain from integration by parts is zero. You can’t say anything about each component of the sum individually.

Nope…! I think you are missing something. If \vec J \cdot \hat n =0 then you cannot conclude that TS\vec J \cdot \hat n =0, and in fact, in my problem, the latter differs from 0. Here is why:

Assume \vec J \cdot \hat n =0 on all boundaries, and in particular on the 2 boundaries with constant r values. On those 2 specific boundaries, \vec J =J_\theta(r, \theta) \hat \theta and \hat n = \hat r, so it is pretty clear that \vec J \cdot \hat n =0 on those 2 boundaries.

Now consider the term S\vec J on those 2 boundaries. S =\begin{bmatrix} S_{xx} & 0 \\ 0 & S_{yy}\end{bmatrix}, or, changing basis to polar, S_\text{polar}= \begin{pmatrix} S_{xx}\cos^2(\theta) + S_{yy}\sin^2(\theta) & -S_{xx}\sin(\theta)\cos(\theta) + S_{yy}\sin(\theta)\cos(\theta) \\ -S_{xx}\sin(\theta)\cos(\theta) + S_{yy}\sin(\theta)\cos(\theta) & S_{xx}\sin^2(\theta) + S_{yy}\cos^2(\theta) \end{pmatrix}.

Now, if you compute S_\text{polar} \begin{pmatrix} 0 \\ J_\theta \end{pmatrix}, you get a vector that has both an \hat r and \hat \theta components. Since T is a scalar field, the flux TS\vec J will possess a non zero normal component (a non zero \hat r component along the 2 boundaries with constant radii). Therefore TS\vec J \cdot \hat n \neq 0 despite the fact that \vec J \cdot \hat n =0.

Again, \vec J_U = -\kappa \nabla T +TS\vec J +V\vec J. Where V\vec J \cdot \hat n = 0 on those 2 boundaries, and \kappa \nabla T may or may not vanish there (it’s up to me, I can choose to apply Dirichlet b.c.s. there or not).

As a conclusion, \vec J_U \cdot \hat n \neq 0 on those 2 boundaries.

Edit: Re-reading the last paragraph of your last post, let me answer to it more clearly. The 2 conditions \vec J \cdot \hat n=0 and \vec J_U \cdot \hat n =0 cannot be fullfilled at the same time on the 2 boundaries where I do not specify any Dirichlet b.c.s. Physically I do know that the first condition is a must, and thankfully FEniCSx correctly applies it in my case. The second condition does not hold, and FEniCSx is able to “see it”, too. Therefore, as is, my code produces a \vec J_U \cdot \hat n \neq 0 on the 2 boundaries where I did not specify any Dirichlet b.c.s.
And about the 2 terms I mention, they are part of \vec J_U, and I focused on TS\vec J in particular so that you can see why \vec J_U \cdot \hat n =0 cannot hold if \vec J \cdot \hat n =0 holds. Hope it’s clearer now.

[Stein]

@raboynics

If \vec J \cdot \hat n =0 then you cannot conclude that TS\vec J \cdot \hat n =0,

Ah right, I overlooked S being a tensor.

@raboynics

The 2 conditions \vec J \cdot \hat n=0 and \vec J_U \cdot \hat n =0 cannot be fullfilled at the same time

I don’t see how that follows from your text.

\vec J_U \cdot n = -\kappa \nabla T \cdot n +T(S\vec J) \cdot n +V\vec J \cdot n = 0
Given J \cdot n = 0, you would get \kappa \nabla T \cdot n = T(S\vec J) \cdot n .

@raboynics

Ok… so, I am sorry but I do not understand how to correctly derive the weak form I seek

Is this problem formulation what you’re trying to achieve?
\begin{cases} \nabla\cdot \vec J_U = 0 \text{ in }\Omega\\ \nabla\cdot \vec J = 0 \text{ in }\Omega\\ J\cdot n = 0 \text{ on }\partial\Omega \\ T = T_{D} \text{ on }\partial\Omega_D \\ -\kappa\nabla T \cdot n= 0 \text{ on }\partial\Omega \setminus \partial\Omega_D \end{cases}
with \vec J_U = -\kappa \nabla T +T(S\vec J) +V\vec J and \vec J =-\sigma \nabla V - \sigma S \nabla T

Hmmm, I’ve puzzled for a bit, but I cannot come up with a clean weak form (without variational crimes) that corresponds to that PDE. To me, that hints at the incompatibility of those BCs, but I could be wrong.

[raboynics]

@Stein

I don’t see how that follows from your text.

\vec J_U \cdot n = -\kappa \nabla T \cdot n +T(S\vec J) \cdot n +V\vec J \cdot n = 0
Given J \cdot n = 0, you would get \kappa \nabla T \cdot n = T(S\vec J) \cdot n .

You’re absolutely right, this time an overlook from my part.

@Stein

Is this problem formulation what you’re trying to achieve?
\begin{cases} \nabla\cdot \vec J_U = 0 \text{ in }\Omega\\ \nabla\cdot \vec J = 0 \text{ in }\Omega\\ J\cdot n = 0 \text{ on }\partial\Omega \\ T = T_{D} \text{ on }\partial\Omega_D \\ -\kappa\nabla T \cdot n= 0 \text{ on }\partial\Omega \setminus \partial\Omega_D \end{cases}
with \vec J_U = -\kappa \nabla T +T(S\vec J) +V\vec J and \vec J =-\sigma \nabla V - \sigma S \nabla T

Hmmm, I’ve puzzled for a bit, but I cannot come up with a clean weak form (without variational crimes) that corresponds to that PDE. To me, that hints at the incompatibility of those BCs, but I could be wrong.

I think you are right about your conclusion.
Ok, thinking more about it, the 4 first conditions are an absolute must. The 5th condition is probably wrong (sorry about that!), I would rather have \vec J_U \cdot \hat n=0 on the boundaries where no Dirichlet conditions are enforced.

I checked up quickly last night (between midnight and 1 am… I was tired to say the least) and I think my current problem is indeed like this, and I was wrong in claiming that \vec J_U \cdot \hat n \neq 0 on the boundaries without any specified boundary conditions (FEniCSx returned numbers like 1e-4 or 1e-5 so I wasn’t sure if this was supposed to vanish or not, but since I used the usual non accurate way to compute fluxes, it may indeed be 0).

Therefore, to sum up, I think I want:
\begin{cases} \nabla\cdot \vec J_U = 0 \text{ in }\Omega\\ \nabla\cdot \vec J = 0 \text{ in }\Omega\\ J\cdot n = 0 \text{ on }\partial\Omega \\ T = T_{D} \text{ on }\partial\Omega_D \\ \vec J_U \cdot n= 0 \text{ on }\partial\Omega \setminus \partial\Omega_D \end{cases} .

This would indeed imply that \kappa \nabla T \cdot \hat n = T(S\vec J) \cdot \hat n on the boundaries where I do not prescribe the Dirichlet b.c.s.

I am getting convinced this is what happens physically. It means that there are 2 boundaries from which energy enters and leave the system, which makes sense since the system is in thermal contact with reservoirs there, and the other 2 boundaries are thermally and electrically insulated. On those insulated boundaries, strangely \nabla T \neq \vec 0, but this is required for no energy exchange with the exterior, because TS\vec J \cdot \hat n \neq 0 there, too, which is a bit weird but apparently possible.

So nice, I think you basically either figured out yourself or at least helped me to understand what is going on in this uncommon system I am dealing with. Things that we take for granted are not necessarily valid anymore, like the thermal insulation implying \nabla T =\vec 0 does not hold in this thermoelectric system.

Thanks a lot Stein, for all the information and your patience.

My question would be, then, whether I am correctly modelling the 5 conditions with my current weak form (with the sigma term corrected)? I think so but I prefer to be 100% sure.

[Stein]

@raboynics

My question would be, then, whether I am correctly modelling the 5 conditions with my current weak form (with the sigma term corrected)? I think so but I prefer to be 100% sure.

Yep, that’s what your weak form is doing

[raboynics]

Perfect, thanks for all Stein
I’ll probably come back with the fluxes at some point, or possible some other questions.