Issue with the accuracy of the computed Flux

variational formulation
1

Hello,

In the code below I tried to create a simple parallel plate capacitor and calculate the charge on both plates.
The problem I faced is that the charge is not the same on both plates which actually should be the case. I tried locally refining the mesh and it did improve the solution but still there was a 10% drift error and even more on a 3D mesh.

Is there a way to increase the accuracy of the solution other than refining the mesh ???

Thank you in advance

mesh = RectangleMesh(Point(0, 0),Point(400, 50), 80, 10)
V = FunctionSpace(mesh, ‘CG’, 1)
tol = 1e-14
class Up(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and x[1] >= box_height - tol

class Down(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and x[1] < tol
and x[0] >= 100 and x[0] <= 300
up = Up()
down = Down()

bcs = [DirichletBC(V, Constant(5.0), up),DirichletBC(V, Constant(0.0), down)]
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(0)
a = dot(grad(u), grad(v))dx
L = fv*dx
u = Function(V)
solve(a == L, u, bcs)

boundaries = MeshFunction(“size_t”, mesh,mesh.topology().dim()-1)
boundaries.set_all(0)
up.mark(boundaries, 1)
down.mark(boundaries, 2)

ds = Measure(“ds”)[boundaries]
n = FacetNormal(mesh)
m1 = dot(n,grad(u))*ds(1)
Q1 = assemble(m1)
—>Q1= 24.75939799680508

m2 = dot(n,grad(u))*ds(2)
Q2 = assemble(m2)
—>Q2 = -22.287241237050743

2

what i mean is that the charge on the upper plate Q1 should be equal to the charge on the lower plate Q2. but this is not what happening

3

The “obvious” approach of directly evaluating \nabla u\cdot\mathbf{n} on the boundary turns out not to be the most accurate way to post-process fluxes from finite element solutions. You can look at my responses here and here for some discussion on extracting boundary fluxes that exactly satisfy the conservation law (including code examples and references).

1 Like
4

Thank you for answering. i checked both responses of yours.
i forgot to mention that, i am still new user to fenics and found your code hard to follow. if you could explain to me how to redfine the flux in my code or which approch should i take that would be helpful

thank you again for your help

5

The theory behind my answers linked above is to try to solve for a flux such that the Neumann problem prescribing that flux would have the same solution as the Dirichlet problem. So, in traditional mathematical notation, you first solve for u^h\in S^h_\text{BC} such that

\int_\Omega \nabla u^h\cdot\nabla v^h\,d\Omega = \int_\Omega fv^h\,d\Omega~~~\forall v^h\in V^h_\text{BC}\text{ ,}

where S^h_\text{BC} is a trial space satisfying the Dirichlet BCs on \Gamma_D\subset\partial\Omega and V^h_\text{BC} is the corresponding test space of functions going to zero on \Gamma_D. Then, instead of directly looking at \nabla u^h\cdot\mathbf{n} on \Gamma_D, you would solve for a function q^h in the \Gamma_D-trace space of V^h, where V^h denotes the discrete space without any boundary conditions and “trace space” means functions restricted to (a subset of) the boundary. The variational problem for q^h is then: Find q^h\in\operatorname{tr}_{\Gamma_D}V^h such that

\int_{\Gamma_D} q^hv^h\,d\partial\Omega = \int_\Omega \nabla u^h\cdot\nabla v^h\,d\Omega - \int_\Omega fv^h\,d\Omega~~~\forall v^h\in\operatorname{tr}_{\Gamma_D}V^h\text{ ,}

where u^h is taken as given (since it was solved for earlier), \operatorname{tr}_{\Gamma_D}V^h denotes the trace space, and you can see from integration by parts that the unknown function q^h plays the role of \nabla u^h\cdot\mathbf{n}. However, q^h converges more rapidly to the flux from the exact solution, and it will satisfy global conservation exactly (on any mesh). This is a fairly old idea, discussed, e.g., here in the 1980s, but remains surprisingly obscure.

The code examples are admittedly messy, since my implementation of the trace space is a bit of a hack. There is probably a much cleaner implementation possible using @cdaversin 's new mixed-dimensional functionality, but I haven’t gotten around to learning the usage of that yet.

2 Likes
6

thanks a lot for the explanation. now i have an idea of how to do it.

you said the function q has the role of ∇u⋅n. does that mean if i want to assemble the flux on the boundary of intrest
i write at the end

assemble(q*ds(boundary-marker).
.

7

Yes, that’s the idea.

8

Hello again,

now i get better flux on the upper boundary(plate) but at the lower plate the flux is zero. They should actually be the same (Q1 = Q2).

Here is the new code :

from fenics import *
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

mesh = RectangleMesh(Point(0, 0),Point(400,50), 80, 10)
n = FacetNormal(mesh)
V = FunctionSpace(mesh,“Lagrange”,1)

tol = 1e-14

class Up(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and x[1] > 50 - tol \

class Down(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and x[1] < tol
and x[0] > 100 and x[0] < 300 and on_boundary

up = Up()
down = Down()

class AntiBdry1(SubDomain):
def inside(self,x,on_boundary):
return x[1] < 50.0 - tol #and x[1]> tol and x[0] <= 100 and x[0] >= 300
class AntiBdry2(SubDomain):
def inside(self,x,on_boundary):
return x[1]> tol
and x[0] < 100 and x[0] > 300

FLUX_BDRY = 1
FLUX_BDRY1 = 3
COMPLEMENT_FLUX_BDRY = 2
boundaryMarkers = MeshFunction(“size_t”,mesh,mesh.topology().dim()-1,
COMPLEMENT_FLUX_BDRY)

up.mark(boundaryMarkers,FLUX_BDRY)
down.mark(boundaryMarkers,FLUX_BDRY1)
ds_marked = ds(subdomain_data=boundaryMarkers)

BCs = [DirichletBC(V, Constant(5.0), up),DirichletBC(V, Constant(0.0), down)]
antiBCs = [DirichletBC(V,Constant(0.0),AntiBdry1()),DirichletBC(V,Constant(0.0),AntiBdry2())]

f = Constant(0)
def R(u,v):
return (inner(grad(u),grad(v)) - f*v)*dx

u = TrialFunction(V)
v = TestFunction(V)
F = R(u,v)
a = lhs(F)
L = rhs(F)

u_h = Function(V)
solve(a==L,u_h,BCs)

qn = TrialFunction(V)
n = FacetNormal(mesh)
consistencyTerm = inner(grad(u_h),n)vds_marked(COMPLEMENT_FLUX_BDRY)
FBdry_up = qnvds_marked(FLUX_BDRY) - R(u_h,v) + consistencyTerm
FBdry_down = qnvds_marked(FLUX_BDRY1) - R(u_h,v) + consistencyTerm

def solveFor_qn_h(FBdry, BCs):
aBdry = lhs(FBdry)
LBdry = rhs(FBdry)
ABdry = assemble(aBdry,keep_diagonal=True)
bBdry = assemble(LBdry)
[BC.apply(ABdry,bBdry) for BC in BCs]
qn_h = Function(V)
solve(ABdry,qn_h.vector(),bBdry)
return qn_h

qn_h_up = solveFor_qn_h(FBdry_up, antiBCs)
qn_h_down = solveFor_qn_h(FBdry_down, antiBCs)

assemble(qn_h_upds_marked(FLUX_BDRY))
assemble(qn_h_downds_marked(FLUX_BDRY1))

9

I think the main issue here is in the definition of the antiBCs (which is inherently confusing due to the “hacky” nature of implementing trace space functions in the volume function space and, again, could probably be cleaned up a lot using the mixed-dimensional functionalitiy). You want to exclude unknowns away from each of the Dirichlet boundaries individually, not both simultaneously, so you would define two “anti-boundary” SubDomains, like

class AntiBdry1(SubDomain):
    def inside(self,x,on_boundary):
        return x[1] < 50.0 - tol
class AntiBdry2(SubDomain):
    def inside(self,x,on_boundary):
        return x[1] > tol or x[0] < 100 or x[0] > 300

(where I’ve fixed the logic in the second one) and corresponding DirichletBCs

antiBC1 = DirichletBC(V,Constant(0.0),AntiBdry1(),"pointwise")
antiBC2 = DirichletBC(V,Constant(0.0),AntiBdry2(),"pointwise")

using the "pointwise" mode to pin down degrees of freedom not associated with boundary facets. Then you would use one at a time when solving for fluxes on the separate Dirichlet boundaries, i.e.,

qn_h_up = solveFor_qn_h(FBdry_up, [antiBC1,])
qn_h_down = solveFor_qn_h(FBdry_down, [antiBC2,])

instead of applying both in each solve. Also, the extra consistencyTerm is also not necessary here, as it is a trick for extracting fluxes from directly adjacent portions of a Dirichlet boundary, so you can simplify the residuals for the q^h problems to

#consistencyTerm = inner(grad(u_h),n)*v*ds_marked(COMPLEMENT_FLUX_BDRY)
FBdry_up = qn*v*ds_marked(FLUX_BDRY) - R(u_h,v) #+ consistencyTerm
FBdry_down = qn*v*ds_marked(FLUX_BDRY1) - R(u_h,v) #+ consistencyTerm

With all these changes, I’m getting fluxes on the two Dirichlet boundaries that cancel out to within very high precision.

P.S. It’s better to format code using backticks (and optionally specifying a language for syntax highlighting), e.g.,

```python
#
# code
#
```
2 Likes