Evaluating an integral yields a wrong result

General
1

I need to evaluate the volume (or surface because my domain is 2D) integral of -T\left(S_{xx}\frac{\partial J_x}{\partial x} + S_{yy}\frac{\partial J_y}{\partial y} \right). Where T is a function of position, \vec J is a position-dependent vector field.
FEniCSx computed correctly \vec J as well as T, very accurately computed, and I was even able to compute another integral involving |\vec J|^2 and got a match up to around 8 significant digits. However, when I try the following code, I do not get the expect result at all, not even close.

assemble_scalar(form(- T * (S_xx * J_vector[0].dx(0) + S_yy * J_vector[1].dx(1)) * dx))

My dx element is the same I use everywhere else without problem and comes from reading a file.msh:

self.mesh, self.cell_markers, self.facet_markers = gmshio.read_from_msh(mesh_file, MPI.COMM_WORLD, gdim=2)
        self.x = SpatialCoordinate(self.mesh)
        self.dx = Measure("dx", domain=self.mesh, subdomain_data=self.cell_markers)

The expected result of the integral is the one I derived by hand/Simpy. Making a MWE is probably not a good idea for this specific question.
If you happen to have an idea on whether my code is equivalent to performing the volume/surface integral on the whole domain, let me know. If I made a mistake, please let me know!
Thanks a lot.

2

How have you defines Sxx and Syy? If those are interpolated on C0 FEM spaces, that could well be the source of the discrepancy

3

They are just constants:

S_xx, S_yy = (1e-5, 1e-3).

4

Hmm, and how is J implemented?
Could you manually check if the integral of dJx/dx is computed correctly?

5

J comes from the solving of a mixed elements problem where V and T are solved for (these are scalar fields), and is constructed from the gradient of V and T.

It looks correct, and the evaluation of the integral of \rho |\vec J|^2 over the same domain yields the expected result simpy provides.

I also checked that div J equals 0, it is worth something e-10 so reasonably close to 0 which looks OK to me. I don’t know of another way to make sure that dJx/dx is computed accurately. If you have any idea, please let me know.

But you too believe that the integral computed is the one I wrote above? If do… damn, what a problem.

6

Indeed I don’t see anything funny. Just to be sure, you’re computing on 1 core right? Else you’d still have to sum over all cores.

You could try to reimplement this to see if the result stays the same:

ddS = ufl.as_tensor([[ S_xx,0 ],[0,S_yy]])
assemble_scalar(form(- T * (ufl.div( ufl.dot(ddS,J_vector) * dx))
1 Like
7

Making an MWE with pre-determined functions T,J_x, J_y,S_xx, S_yy (through interpolation of spatially dependent functions), would point you to if there is something fundamentally wrong with the formulation above.

It would also greatly simplify the amount of guess-work anyone interested in reproducing this behavior would have to do

1 Like
8

The problem is that there is no known analytical solution for simple geometries, I could possibly derive some, like I have done for a quarter of an annulus, but it might take me months to do so (like it did). My MWE would be gigantic if I provided my code, even trying to reduce it as much as I can. I have written it in OO form, and I read a mesh file from gmsh (quarter of an annulus I cannot simplify for the MWE because I would lose the analytical formula).

I don’t mind doing it, let me know if you believe it might help.

9

Yes 1 core. I will try your suggestion tonight, thanks.

10

It wouldn’t have to be an analytical solution to your problem. @dokken is suggesting you simply define some vector valued J and some scalarvalue T for which you can work out that integral, and then see if it matches FEniCS’ output.

1 Like
11

Ah, I think I see what you mean. Yes, I can get the analytical solution for T, Jx and Jy, even their derivative. And I could compare with the one computed by,fenicsx, I am just not fluent enough with fenicsx to do so.

I guess I have to define J_analytical, interpolate in some high degree space, and compute the derivatives… and then visually inspect the difference between the computed Jx and the analytical one.

Is this what you have in mind?

12

Yes, that is what I had in mind.

Please note that the variational forms does not differ much between DOLFINx and legacy FEniCS, and you should be able to verify that the legacy version of an analytical example gives the same as the DOLFINx one.

1 Like
13

Thanks a lot guys. It turns out I had made at least a mistake in the analytical formula, correcting it yields an expected result of 0.012850 whereas FEniCSx yields 0.012780. It isn’t nearly as good as the other integral (I get like 8 significant digits matching), but I guess it’s fine.

I do not run the legacy version of FEniCS anymore, I’d rather not convert my code to legacy fenics :slight_smile: due to lack of time.

I still get non sensical results for some quantites such as div J_U integrated over the whole volume, which should yield 0 but yields -2.87. But it’s another topic I suppose. In any case I should learn how to verify the numerical solution with a given analytical expression… I will try when I get the time.
Thanks again guys. I will probably come back to you after my attempt.