Computing a flux involving 2 functions (ME)

General
1

As a follow up to Computing fluxes accurately, where I created a MWE to compute a flux and eventually reached a nice result using the method described by Jeremy Bleyer (Computing consistent reaction forces — Computational Mechanics Numerical Tours with FEniCSx), I face several problems or questions when I apply this method in my real problem.
It is a mixed element problem where I solve for two scalar fields, T and V.
I build up a vector field \vec J=-\sigma \nabla V -\sigma S \nabla T that involves the 2 unknown scalar fields, and I write the weak form of 2 coupled PDEs, the PDEs being \nabla \cdot \vec J =0 and \nabla \cdot (\kappa \nabla T)+\vec J^2/\sigma-T(S_{xx}-S_{yy})\frac{\partial J_x}{\partial x}=0.

I am able to solve for V and T, but then I need to compute fluxes, several kinds of fluxes, on each boundary of the domain. And this is where I have several problems. For example, computing the heat flux through each boundary does not yield very accurate results, where I expect 0, I can get things that differ significantly from 0. And it’s not just fluxes, but also when I integrate the divergence of those fluxes in the whole volume, I do not exactly recover Gauss’s law, i.e. a volume integral of the divergence of a field can be significantly different than the surface integrated flux of that field, which is very bad.

I have several questions.
If I use Jeremy’s method on a boundary, which flux is computed exactly in my problem? In my MWE, it was obvious it was the heat flux, but here, I have a ME space, I have no idea which flux is being computed with this code:

(...)
# Define ME function space
deg = 4
el = FiniteElement("CG", mesh.ufl_cell(), deg)
mel = MixedElement([el, el])
ME = FunctionSpace(mesh, mel)
u, v = TestFunction(ME)
TempVolt = Function(ME)
temp, volt = split(TempVolt)

# Weak form.
J_vector = -sigma * grad(volt) - sigma * Seebeck_tensor * grad(temp)

F_T = dot(-κ * grad(temp), grad(u)) * dx + dot(rho * J_vector, J_vector) * u * dx - temp * (s_xx * J_vector[0].dx(0) + s_yy * J_vector[1].dx(1)) * u * dx

F_V = -dot(grad(v), sigma * grad(volt))*dx -dot(grad(v), sigma * Seebeck_tensor * grad(temp))*dx - v * J_out * ds(out_current_curve) + v * J_inj * ds(inj_current_curve) 

weak_form = F_T + F_V
(...)

# Compute heat flux on boundaries more accurately than usual method.
u_bc = fem.Function(ME)
bcs_right = [fem.dirichletbc(u_bc, right_dofs_temp)]    
virtual_temp_right = Function(ME)
u_bc.sub(0).interpolate(one)
fem.set_bc(virtual_temp_right.vector, bcs_right)    
virtual_heat_flux_form_right = form(action(weak_form, virtual_temp_right))
virtual_heat_flux_right = assemble_scalar(virtual_heat_flux_form_right)     
print(f'Virtual flux right computation: {virtual_flux_right}')

The above code returns a number very close to what I obtain by the standard way to compute fluxes, except that I do not know which flux is computed.

Since the heat equation I am solving is in fact worth \nabla \cdot \vec J_U=0 where \vec J_U=-\kappa \nabla T +(ST+V)\vec J is the internal energy flux, I suspect Jeremy’s method gives me the internal energy flux over the boundary I specify. But it’s not clear at all, to me. The “sub(0)” part annoys me somehow.

What happens for example if I do:

u_bc.sub(1).interpolate(one)

either in addition to u_bc.sub(0).interpolate(one) or on its own?

Ideally, I wish to compute with utmost accuracy the \vec J_U flux across each boundary. And if possible, also \vec J. The usual method is not accurate enough, even when I use a space of degree 7 (Jeremy’s method converges much better than the usual method, but I do not know how to compute a particular flux with it).

My MWE would be around 300 lines of code. I don’t mind sharing the code, but maybe Bitbucket would be more appropriate, if there is any need to share the full code…

2

If weak_form is a linear form L(w), this line essentially computes L for a given w. In your case, w = (u,v) and L(w) = F_T(u)+F_V(v). Then, you choose w=u_{bc} and set the sub(0) component to the function one, meaning w = (one,0)). So, this is computing L( (one,0) ) = F_T(one) + F_V(0) = F_T(one).

Then you would compute L( (one,one) ) = F_T(one) + F_V(one) or L( (0,one) ) = F_T(0) + F_V(one) = F_V(one) , respectively.

Does that help you figure out your remaining questions?

2 Likes
3

This helps a lot, yes.

I do not understand my results, however. For example with u_bc.sub(0).interpolate(one), I get a value very close to the heat flux I can compute, which happens to be also very close to the energy flux too, so I cannot really distinguish between which one is reached. (the values are about 740 if I remember well).

When I do u_bc.sub(1).interpolate(one) however, I reach “0.0”, not 1e-9 or something else that is close to 0. I get exactly 0.

The F_T part is the heat equation (which should be better renamed as internal energy flux equation since it is worth \nabla \cdot \vec J_U =0 which is NOT equal to the divergence of the heat flux in general).

The F_V part is \nabla\cdot \vec J =0, and I still get 0.0 even when I impose a flux with Neumann boundary conditions of I = 2 A (the current equals 2 amperes), so I know I should get exactly 2.0, but I get 0.0.

I tested many variants… never could get anything else than 0 for the J flux.

Therefore I do not even understand what I get as flux when I use the sub(1) part…

4

It’s really just a matter of working out the linear form to see what you get:

F_T(u) = \int_\Omega \Big[ -κ \nabla T \cdot \nabla u + \rho J \cdot J u - T ( s : \nabla J ) u \Big] dx
F_V(v) = \int_\Omega \Big[ - \sigma \nabla V \cdot \nabla v - (\sigma S \cdot \nabla T ) \cdot \nabla v \Big] dx - \int_{S_{out}} v J_{out} ds + \int_{S_{inj}} v J_{inj} * ds

Looking more closely, I now understand that you’re constructing virtual_temp_right as the function that is zero everywhere and interpolates 1 at the dirichlet boundary. Plugging in that function for u (let’s call it q) and performing integration by parts, gives:

\int_\Omega \Big[ \nabla\cdot(κ \nabla T )+ \rho | J | ^2 - T ( s \cdot \nabla J ) \Big] q \, dx - \int_{S_{dirichlet}} κ \nabla T\cdot n \, ds
where the \nabla\cdot(κ \nabla T ) term must be understood distributionally. The first integral is the residual of de PDE, limiting to zero. The second term is the flux extraction.

In your second case, I wonder if the Dirichlet boundary for T is also the Dirichlet boundary for V? Else, v=q is simply in the test-space and by definition F_V(q) = 0.

P.S. Your term \int T ( s : \nabla J ) u \, dx is formally not permitted, since it effectively makes use of second derivatives of functions that are only first order continuous. That’s called a variational crime and could cause issues, especially if you’re so focused on accuracy and convergency. You’d have to perform integration by parts somehow.

1 Like
5

Thanks a lot Stein, I will try to fully understand everything you wrote.

What I can say is that, as of now, no, there is no Dirichlet b.c. for “volt”, or V. I let the solver determine the “zero” of V, since it is arbitrary. I did not face any problem so far, the code converges.

I am a bit worried about what you said about the second derivatives. Yes I realize that I compute \frac{\partial J_x}{\partial x} where J_x involves both \frac{\partial T}{\partial x} and \frac{\partial V}{\partial x}, and so there is indeed a second order derivative. Anyway… this gets me worried a lot. I do not get any error, and the results matches closely an analytical computation of a similar but simplified problem.

However I remember dokken mentioning that my \vec J is actually discontinuous even though V and T are continuous. So, are you sure that I am dealing with functions that are first order continuous? I would have thought 0th order continuous.

If I understand you well, my result computes \int _{S_\text{Dirichlet}}\kappa \nabla T \cdot n ds when I use sub(0).

Maybe I could define smartly a new residual form in order to extract \int \vec J \cdot \hat n ds over any boundary, and the same for \vec J_U?

I thank you again Stein for all of this… I am not in the world of finite elements, there is a lot to learn for me.

6

Then indeed any function you can choose for v by definition satisfies F_V(v)=0. That’s what the discrete formulation enforces.

T and V are C0 continuous. That means J \approx \nabla T is discontinuous. That means \nabla J is taking a derivative of a discontinuous function. That is a variational crime.

Yes, that’s correct.

Typically the fluxes are known at all boundaries except for the Dirichlet boundary though? For J, on S_{out} you prescribe J_{out}, on S_{inj} you prescribe J_{inj}, and on all the rest you prescribe 0. For J_U = -\kappa\nabla T you prescribe 0 on all but the Dirichlet boundaries.

1 Like
7

Thanks once again Stein for your informative reply!

Not really in my case. I am dealing with a thermoelectric problem. I can use Dirichlet boundary conditions on 2 boundaries for the temperature variable (it’s a 2D geometry), and this settles every fluxes. I can also enforce an entering \vec J flux through any boundary I wish, and if I do that, yes the total net flux of \vec J through each boundary will be known to me. But not \vec J_U, which by the way is worth -\kappa \nabla T + (ST + V)\vec J and not just -\kappa \nabla T.

By imposing a given temperature on 2 boundaries, everything adapts/settles, the problem is well posed (with the caveat of the voltage’s zero). The fluxes need to be computed though, in general. In my case they are too hard to compute exactly by hand, and I actually use FEniCSx to guide me to help me understand the physics of the problem.

I’ll now attempt to solve correctly my problem and deal with the -T(S_{xx}\frac{\partial J_x}{\partial x} + S_{yy}\frac{\partial J_y}{\partial y}) term using integration by parts.

8

I’ll have to double check the computations.

The term -T(S_{xx}\frac{\partial J_x}{\partial x} + S_{yy}\frac{\partial J_y}{\partial y}) is worth -T(S_{xx}-S_{yy})\frac{\partial J_x}{\partial x} because \nabla \cdot \vec J=0 holds.

Using your notation of the test function q instead of u for the heat equation, integrating the above term multiplied by q, using integration by parts, I reach -\sigma\int_\Omega \left( \frac{\partial V}{\partial x} + S_{xx}\frac{\partial T}{\partial x} \right) (S_{xx}-S_{yy}) (q \nabla T +T \nabla q)d\Omega.

Therefore I would replace my code - temp * (s_xx * J_vector[0].dx(0) + s_yy * J_vector[1].dx(1)) * u * dx which is the last term of F_T(q) by the expression I just derived. I do hope this problem is well set now, assuming I didn’t make any mistake?

I assumed that the surface/boundary term involving q vanishes.

9

Hmm unfortunately, I get an error when I use this code:

- (( sigma* volt.dx(0) + sigma * s_xx * temp.dx(0))*(s_xx - s_yy) * (u * grad(temp) + temp*grad(u))) * dx

The error is

Traceback (most recent call last):
  File "/root/shared/debug_quarter_annulus/test.py", line 528, in <module>
    post_proc = PostProcessor(thermoelec_sim)#thermoelec_sim, display_output=True, save_J_e=True)
  File "/root/shared/debug_quarter_annulus/test.py", line 239, in __init__
    self.thermoelec_sim_res = thermoelec_sim.run_sim()
  File "/root/shared/debug_quarter_annulus/test.py", line 179, in run_sim
    newly_derived_term = - (( sigma* volt.dx(0) + sigma * s_xx * temp.dx(0))*(s_xx - s_yy) * (u * grad(temp) + temp*grad(u))) * dx
  File "/usr/local/lib/python3.10/dist-packages/ufl/measure.py", line 361, in __rmul__
    raise ValueError(
ValueError: Can only integrate scalar expressions. The integrand is a tensor expression with value shape (2,) and free indices with labels ().

I do understand that I am apparently trying to integrate a vector, which shouldn’t be the case… hmm. I guess I did a mistake in the weak form…

Edit: I cannot find any mistake. Not sure why I can’t implement the weak formulation.

Edit 2: I could bypass the problem by rewriting the full weak form for F_T as follows: dot(-κ * grad(temp), grad(u)) * dx + dot(sigma * temp * Seebeck_tensor* J_vector, grad(u)) * dx + dot(volt * J_vector, grad(u)) * dx. I hope that under this form I do not have 2nd order spatial derivatives…
It yields almost exactly the same results as before, roughly about 1e-7 difference than the weak form I started this thread with.
However I think this changed a bit the fluxes computed using the accurate way. I would have to figure this out tomorrow.
Thanks for all, and if you spot anything wrong please let me know.

10

I must admit that I am getting a little lost in the different flux definitions, so I won’t dare to make any definitive claims, but looking at your weak form it seems to me you need two BCs at every point in the boundary. Either a Dirichlet condition on T and/or V, or a natural condition on J and/or \kappa\nabla T \cdot n .

If you don’t add an integral term in place of J or \kappa\nabla T \cdot n on the boundary, but you also don’t specify a Dirichlet condition at that boundary, then you are implicitly saying J=0 or \kappa\nabla T \cdot n=0.

The form looks valid to me (I haven’t checked your derivation steps though).

That’s great! Indeed, the variational crime could have been small, so one would not expect big differences, but it is always better to do things properly.

Note that now you have changed your weak form, and exectured more integration by parts steps, you have to revisit the meaning of your flux extraction…

1 Like
11

Hello Stein and thanks for replying.

I might create a new thread just to make sure my weak form corresponds to the PDEs I have in mind, I have a little doubt about the flux terms. Because I just remembered that the test function vanishes only where Dirichlet b.c.s are applied, and I have 2 boundaries without any Dirichlet b.c.s, therefore I wonder if I am missing a flux term.

But… I am sure that prescribing only 2 temperatures on two boundaries is enough to fully determine the system. The temperature difference creates a current density inside the material, \vec J=-\sigma \nabla V - \sigma S \nabla T. Where V adjusts itself on all boundaries so that \vec J is parallel to all boundaries. So, if I just apply 2 Dirichlet b.c.s., this is what I expect, and this is also what I get with FEniCSx (and what I got analytically in a simplified problem where there is no intermixing of V and T).

The “sad” part is that no flux is trivial, IMO, so I need FEniCSx to guide me, and I need a very accurate flux evaluation because, unfortunately, integrating \vec J_U along a boundary yields very similar results than integrating -\kappa \nabla T there, and I need to have an accurate evaluation of their differences.

I think I will reach it, I will study again how you derived that I had computed -\kappa \nabla T flux. To be honest, I couldn’t follow the first part, where you apparently replaced u by 1 and integrated by parts to reach \int_\Omega \Big[ \nabla\cdot(κ \nabla T )+ \rho | J | ^2 - T ( s \cdot \nabla J ) \Big] q \, dx - \int_{S_{dirichlet}} κ \nabla T\cdot n \, ds. Why is the surface integral only involving \kappa \nabla T?

12

On those boundaries, you’re implicitly prescribing that the (normal components of the) fluxes are zero.

Getting a well-posed system is 1 thing. Whether it relfects the reality that you’re trying to mirror is another :wink:

That’s the term you get from reverse integrating by parts the -\int_\Omega \kappa\nabla T\cdot \nabla q term in your weak form.

1 Like
13

Hmm, which fluxes exactly? I do wish \vec J to have no normal component, yes. But there are many other possible fluxes involved, for example \vec J_U. I expect this one to have a non zero normal component, at least in one configuration of my problem.

The region is 2D, a quarter of an annulus where \theta ranges from 0 to \pi/2, and if I apply Dirichlet boundary conditions for the temperature at the \theta = 0 and \theta =\pi/2 boundaries, I do expect a non zero \vec J_U entering/leaving at the other 2 boundaries where no Dirichlet b.c.s are applied.
You can see that this holds because \vec J_U = -\kappa \nabla T +ST\vec J + V\vec J, with \vec J being parallel to all boundaries. Since V is a scalar field, the term V \vec J contributes nothing, i.e. no entering/leaving flux due to this term. The term ST\vec J however is very different, S being a tensor represented in matrix form as \begin{bmatrix} S_{xx} & 0 \\ 0 & S_{yy}\end{bmatrix}, if you convert it under polar basis, the matrix isn’t diagonal anymore and S\vec J (where \vec J along those 2 boundaries has no radial components, so is worth J_\theta \hat \theta) contains both a \hat r and \hat \theta components, along the 2 boundaries with constant r, meaning that there is a local energy flux that either enters or leaves the material through those 2 boundaries that lack any Dirichlet condition. And there is even a net flux, even when integrated, over those boundaries…

This is why I am confused about how I specify my problem with FEniCSx. Looking at my weak form for example, it is not clear to me at all which condition I apply. Am I applying \nabla T \cdot \hat n= 0 on the boundaries with no Dirichlet boundary conditions? Am I applying, also, \nabla V \cdot \hat n=0 on those boundaries as well?

I guess not.

By the way, thanks a lot for all the rest of your comments I do/did not reply too, they are informative as well, very appreciated.

14

Hello Stein and anyone else reading/interested,
Just so you know, I have created a thread about my weak form (Deriving the weak form of my problem, is it correct?). I think I have understood a bit more what happens when I don’t specify any boundary condition on some boundaries, and to answer my own question, what happens isn’t that any flux vanishes there, but that \nabla T and \nabla V have no normal component there. But a flux, for example TS\vec J may have a normal component there.

I am 99% sure my weak form is correct, and that the results are making a lot of sense.

Now regarding

Indeed, I will try to figure this out myself. I wonder also, is there a smart way to retrieve any flux using the accurate way, not just the one that pops up from the weak form? If I had to choose right now, I’d be interested in \vec J_U, \vec J and TS\vec J as well. Should I try to figure out what kind of weak form yields those fluxes?