Help with discretizing the weak form in time

variational formulation
1

Hi,

I want to solve the time dependent problem of current conduction through media with both conductivity and polarization using FenicsX (0.4.1, conda). The equation is written as

image

After introducing the unknown scalar potential, the final strong form becomes

image

Solving the static solution gives perfectly reasonable results. The problem arises when I try to solve it in time and I would appreciate any help in finding my mistake.

Through backward Euler discretization of the time derivative, we get the following bilinear and linear forms

image

image

In Fenics, I implement it as:

Omega = dolfinx.fem.FunctionSpace(mesh, ("CG", 1))
V = ufl.TrialFunction(Omega) # solution
Vt = ufl.TestFunction(Omega)

U = dolfinx.fem.Function(Omega) # U is the unknown solution
U0 = dolfinx.fem.Function(Omega) # existing/initial solution in time stepping
U.x.array[:] = 0 # initial value
U0.x.array[:] = 0 # initial value

# sigma, epsilonr are domain dependent constants; dt is time step and eps0 is a scalar constant
r = ufl.SpatialCoordinate(mesh)[0] # radial axis
aL = (ufl.dot(ufl.grad(Vt), sigma*ufl.grad(V))*dt + ufl.dot(ufl.grad(Vt), eps0*epsilonr*ufl.grad(V))) * r * ufl.dx
L = ufl.dot(ufl.grad(Vt), eps0*epsilonr*ufl.grad(U0)) * r * ufl.dx

Solving this problem in Fenics exactly following this tutorial gives me a strange result for the potential difference between two points:

image

While the answer does not match with analytical results, the first anomaly I am struggling with is the sudden jump in potential that always happen at the first time step I solve on. I have tried playing with the time steps/interval but this jump is permanent. I implemented the weak form in Comsol (which does time discretization automatically under the hood) and it works well. So the problem is most likely when I do time discretization myself in the variational form. Before I prepare and share an MWE, I wanted to check if people could spot a mistake with the weak form I have?

Whoever wants to help, let me know if I should add more details to make it comprehensible. Thanks a lot.

2

I am now suspecting the problem lies with the boundary conditions. I have direchtlet boundaries at the top and bottom surface of the rectangular field of view. The initial condition of the unknown potential defined to be zero is naturally not compatible with this choice. Does anyone have suggestions on how these cases are solved? I tried to project the initial solution to l2 space as suggested in the tutorial on heat diffusion. However, that only removes the jump by finding an initial state that removes the discontinuity at t=0. Any pointer will be highly appreciated.