Finding the weak form of an eq. with DG elements (finite elements)

[fluidistic]

My goal is to numerically solve (finite elements using FEniCSx software) Ohm's law $\vec J = \sigma \vec E$, where $\vec E = -\nabla V$, $\vec J$ is given (the current density is given on some boundaries), and $\sigma$ is algo given (the electrical conductivity). The problem is solved for $V$, the electrostatics potential, and this allows to compute the electric field in all the space. I can already do that in some region/geometry in space. Now I am attempting to consider 2 stacked materials, in contact with each other, and I want to simulate an electrical resistance between them, by imposing a discontinuity in the electrostatic potential at the (internal) surface where the materials meet.

And this is where I'm having troubles. The weak form of the problem when dealing with discontinuous Galerkin elements (the usual weak form can be found by taking the divergence of the PDE, multiplying it by a test function $v$ and integrating it over the whole domain, applying integration by parts) is apparently very complicated compared to the case where the elements are continuous. For example, the lines 86 to 93 of https://bitbucket.org/fenics-projec...mo/undocumented/dg-poisson/demo_dg-poisson.py show such an example. There are terms from Nitsche that apply to the boundaries where Dirichlet conditions are present, and other terms applied on the interface between the materials where the discontinuity is imposed. But I do not quite understand all of it, and I am not sure how to apply this to my equation.

I have tried to follow the tutorial at https://github.com/jorgensd/dolfinx-tutorial/blob/dokken/dg_tutorial/chapter1/dg.ipynb, but I am still unsure on how to apply it to my equation. If someone could guide me on how to find the expression for the weak form, I'd be happy.

 

[Paul Colby]

Wow, FEniCS is quite neat. Looks like what you need is a variational statement of your problem. It’s a static problem so both $\dot{B}$ and $\dot{D}$ are zero in Maxwells equations. These give us $$E = -\nabla V$$ and $$\nabla \cdot J=0$$. Ohms law gives us, $$E = \sigma J$$ making our equation, $$\nabla\cdot\left(\sigma \nabla V\right)=0$$ The functional you seek is, $$L=\int \sigma |\nabla V|^2 dx^3 $$

 

[fluidistic]

Wow, FEniCS is quite neat. Looks like what you need is a variational statement of your problem. It’s a static problem so both $\dot{B}$ and $\dot{D}$ are zero in Maxwells equations. These give us $$E = -\nabla V$$ and $$\nabla \cdot J=0$$. Ohms law gives us, $$E = \sigma J$$ making our equation, $$\nabla\cdot\left(\sigma \nabla V\right)=0$$ The functional you seek is, $$L=\int \sigma |\nabla V|^2 dx^3 $$

Yeah I know about this.
It's a little bit more involved, the usual weak form involves a test function, denoted by "v" sometimes in the FEniCS manual. I am able to retrieve the weak form for the usual continuous Galerkin elements. I have no problem with this. I can already solve several kinds of problems with FEniCS. E.g. when the materials are anisotropic for example (so sigma is a tensor in this case, rather than a scalar), when there are several materials involved, etc.

I am just stuck when it comes to discontinuous Galerkin elements, in view of imposing a potential drop across the interface of 2 different materials. If you check the links I provided, you would see that the weak form looks like so:

a = dot(grad(v), grad(u))*dx \
   - dot(avg(grad(v)), jump(u, n))*dS \
   - dot(jump(v, n), avg(grad(u)))*dS \
   + alpha/h_avg*dot(jump(v, n), jump(u, n))*dS \
   - dot(grad(v), u*n)*ds(1) \
   - dot(v*n, grad(u))*ds(1) \
   + (gamma/h)*v*u*ds(1)
L = v*f*dx - u0*dot(grad(v), n)*ds(1) + (gamma/h)*u0*v*ds(1) + g*v*ds(2)

for an equation that's not that much different than mine. It is much more complicated than the case where continuous elements are used (and where I have no problem with).

Just to be clearer, I am stuck at finding the weak form that should resemble a lot to the one I post here. It's not clear at all to me where many terms come from.

 

[Paul Colby]

My understanding the weak form of a differential equation is obtained by multiplying the equation through by a set of testing functions then integrating this product over each element. The Galerkin method is to choose the testing functions and expansion functions as the same set of functions. The functional I gave will yield the Galerkin system matrix less the boundary conditions (AKA the right hand side vector).

In this functional, $\sigma$, the conductivity of the materials, is specified. No derivatives of $\sigma$ appear in the functional so there are no restrictions on continuity. Also required is that the potential, $V$, be continuous and differentiable throughout the region. For any real materials, this is exactly what Maxwell's equations dictate. Failing to do this would spawn infinite $E$-fields and we wouldn't want that to happen. So, I don't follow the physics of what you intend interior to the region.

Now, I didn't give the part of the functional that yields the boundary conditions. I believe these are what account for the extra terms seen in the python code. They appear to be sums over the boundary and likely yield a properly tested boundary value equations. For specifying the value of $V$ on the boundary things are conceptually simpler. How this is done in FEniCSx?

 

[Paul Colby]

Okay, I think I see how to do this. Let $R$ be the region and $\partial R$ its boundary. Clearly, our continuum equation is, $$\nabla\cdot(\sigma \nabla V) = 0.$$ We multiply by a testing function $v$ and integrate by parts, $$\int_R \nabla \cdot (\sigma\nabla V)dx^3 = \int_{\partial R} v (J_o\cdot \hat{n})dx^2 - \int_R \sigma \nabla v \cdot \nabla V dx^3=0$$ So, the surface integral is your weakly tested boundary condition where $J_o$ is the specified current on the boundary and $\hat{n}$ the outward directed normal of the boundary.

 

[fluidistic]

My understanding the weak form of a differential equation is obtained by multiplying the equation through by a set of testing functions then integrating this product over each element. The Galerkin method is to choose the testing functions and expansion functions as the same set of functions. The functional I gave will yield the Galerkin system matrix less the boundary conditions (AKA the right hand side vector).

In this functional, $\sigma$, the conductivity of the materials, is specified. No derivatives of $\sigma$ appear in the functional so there are no restrictions on continuity. Also required is that the potential, $V$, be continuous and differentiable throughout the region. For any real materials, this is exactly what Maxwell's equations dictate. Failing to do this would spawn infinite $E$-fields and we wouldn't want that to happen. So, I don't follow the physics of what you intend interior to the region.

Now, I didn't give the part of the functional that yields the boundary conditions. I believe these are what account for the extra terms seen in the python code. They appear to be sums over the boundary and likely yield a properly tested boundary value equations. For specifying the value of $V$ on the boundary things are conceptually simpler. How this is done in FEniCSx?

I am sorry, I am in a rush right now to give you a more elaborated reply. V needs not be continuous. If you place 2 materials next to each other, I can assure you that there will be an electrical resistance across the interface (as well as a thermal resistance, by the way, so the temperature profile isn't continuous either). There is nothing shocking if E is infinite with Maxwell equations, take the E field of a point charge where the charge lies, or the one of a charged infinitely thin wire, right at the wire.

Again, I have no problem deriving.the correct weak form of the eq. I deal with in the case of continuous elements (i.e. what you did, and from memory I think you either got it or got pretty close). But that's not what I wish to achieve.

 

[Paul Colby]

I’ve found this wiki helpful,

Discrete Galerkin Method

Likely material you know. The time dependent current flow is discussed. Superficially this looks similar to the continuous development I sketched, but there is a significant modeling difference. Whereas I assumed a continuous potential function, $V$, this is not strictly necessary. One may use discontinuous vector functions to model $J$ - provided one enforces charge or current conservation at each volume element. This is done by summing the charge flow across element boundaries.

I believe that the terms you seek are just this, the component of current normal to material boundaries must be conserved.

So, this leaves me with the same problem you are faced with. Just replacing $V$ and its derivatives with $J$ I lose the material properties $\sigma$. Varying conductivity must alter the current conservation equations at element boundaries