Comment my code + many questions

[dokken]

@raboynics

One of my goals now is to use the code on different meshes for example, without tweaking the code. What would be a nice way to accomplish this task? Make a function out of my whole code, where as input I would manually specify the physical boundaries where boundary conditions are applied? It would become pretty complex if it could deal with 1D, 2D and 3D meshes without crashing.

I would suggest using class structures.
Then you can have different functions applying appropriate boundary conditions.
Also keep in mind that the only thing that should differ in the core of your variational formulation is the mesh/function space, the UFL forms should be the same in 1D/2D/3D.

Oasisx (disclaimer Im the author) GitHub - ComputationalPhysiology/oasisx: Next generation Open Source Navier Stokes solver using FEniCSx
is an example of how to re-use the variational formulation for different problems.

@raboynics

• As it is currently, Newton’s method fails to converge if I modify Seebeck coefficient to values higher than about 30e-6. However in real materials this number can reach 1000e-6, so I would like to see what would happen in that case, but I am unable to make the Newton’s method converge. Any pointer is welcome.

See for instance: Default absolute tolerance and relative tolerance - #4 by nate

@raboynics

I wanted to specify an electric current magnitude (say 1 ampere), which I used to do by applying a 0 voltage Dirichlet boundary condition on one surface (where the current is injected), and by adding a term in the weak form for the surface where the current exits the material (Neumann’s boundary condition). It used to work well, however it doesn’t seem to work here. I use mixed elements, so I suspect this is the cause, but still, this sounds strange to me. To fix the problem, I considered pure Neumann’s boundary conditions, so I don’t use any Dirichlet boundary conditions for the voltage. As a result, the code seems to find some “zero” of voltage by itself. Voltage is usually defined up to a constant. However the solution I get if I set a boundary to a constant is wrong. But when I let FEniCSx decide which constant it is, it works fine. I do not understand why it is the case. Any pointer is welcome. It is as if FEniCSx thinks the voltage isn’t defined up to a constant.

I would strongly suggest you try to simplify your problem to illustrate this issue.
If you won’t simplify it, you should really be careful about the integration by parts of your problem, making sure that you are not missing a sign in the transformations. Additionally, when you say the code seem to find some “zero” voltage by itself, you should look at the equation you are now solving.

@raboynics

would like to reuse the code with some modularity. For example, I would like to have the choice to ignore thermoelectricity (so by setting Seebeck coefficient to 0). I also want to be able to ignore Joule heating. I can easily do this if I define a function of my whole code and where I have a boolean variable for Joule_heat and where if Joule_heat ==True, then I would include this term in the weak form, otherwise set it to 0 or do nothing. But, as a non programmer, would this be a nice way to achieve this? I realize I do not use any class nor dictionary in my whole code. I am surely missing improvements.

Again, using classes, and boolean inputs to either:

1. Add terms to your variational form
2. Add a c=dolfinx.fem.Constant(mesh, value) to your variational form (i.e. F += c*joule_heating_terms), where value is 0 or 1 depending on if you want to add the effect or not.

I do not have time to go through your code in detail, but I would strongly suggest you to start small, and take out simple parts of your code that you can re-use, and make up a general structure as you go.