Degrees and continuities in FE spaces of several linked quantities

General
1

Hello and sorry for a badly formatted text (cell phone).

In brief, I have a code that solves for a temperature and a voltage (mixed elements, and can be found in my post history, but I’ll probably make the link once I am in front of my computer).

From ‘‘temp’’ and ‘‘volt’’, I compute J and Jq, the current density and the heat flux. They involve taking.the gradient of temp and volt. J appears in the heat eq. as a source term.

From J, I compute dJx/dx and similarly for the ‘‘y’’ component, I add those terms to form a source of heat in the heat equation.

Here are my questions:

Physically, the source term is non zero. However, I get that it vanishes if temp and volt are defined in CG spaces of degree 1, J and Jq are then defined in DG spaces of degree 0, and since I then take the spatial derivatives of a cell wise constant field, I get 0. This makes some sense.

Therefore, if I,want to compute the source term accurately, I need to increase the degrees of all FE spaces. And I don’t know if I can choose a CG space of degree 2 for temp and volt. Because then I,could take CG spaces of deg 1 for J and.Jq, but then this would make a space DG of degree 0 for the source term. And this source term is used to compute temp and volt. In other words, a discontinuous Galerkin space would be used to get a solution that would be continuous Galerkin. Is that possible? Is this mathematically allowed?

If not, does this mean that if I want a real estimate or computation of the solution to the coupled PDEs, I should use DG for all spaces?

2

Here is some relevant code, and other questions.


# Define ME function space
deg = 2
el = FiniteElement("CG", mesh.ufl_cell(), deg)
mel = MixedElement([el, el])
ME = FunctionSpace(mesh, mel)

u, v = TestFunction(ME)
TempVolt = Function(ME)
temp, volt = split(TempVolt)
dTV = TrialFunction(ME)

rho = 1
sigma = 1.0 / rho
κ = 1.8

# The Seebeck coefficient has to be an anisotropic tensor for the Bridgman heat to take place.
s_xx, s_yy = seebeck_components
Seebeck_tensor = as_tensor([[s_xx, 0], [0, s_yy]])

# Skipping the boundary conditions to save space here.
# Weak form.
J_vector = -sigma * grad(volt) - sigma * Seebeck_tensor * grad(temp)
Fourier_term = dot(-κ * grad(temp), grad(u)) * dx
Joule_term = dot(rho * J_vector, J_vector) * u * dx
Bridgman_term = - temp * (s_xx * J_vector[0].dx(0) + s_yy * J_vector[1].dx(1)) * u * dx
F_T = Fourier_term + Joule_term + Bridgman_term
F_V = -dot(grad(v), sigma * grad(volt))*dx -dot(grad(v), sigma * Seebeck_tensor * grad(temp))*dx - v * J * ds(out_current_curve) + v * J * ds(inj_current_curve) 

weak_form = F_T + F_V

# Solve the PDEs.
n, converged = solver.solve(TempVolt)

As you can see, I have a term involving J_vector[0].dx(0) where J_vector involves grad(volt) and grad(temp). Therefore, since I picked a CG degree 2 FE space for temp and volt, I suppose this term in the weak form, internally in FEniCSx, is DG of degree 0. Can someone confirm this?

My other concern is that I need to compute some quantities from the solution. These quantities involve precisely this term of the weak form, integrated over the whole volume. Here is exactly how I do it:
heat_source_term = assemble_scalar(form(- temp * (s_xx * J_vector[0].dx(0) + s_yy * J_vector[1].dx(1)) * dx(domain=mesh))). I need to compute this quantity as precisely as possible. So far, I am not successful in doing so.

As a side note, my solver options are:

problem = NonlinearProblem(weak_form, TempVolt, bcs=bcs, J = Jac)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-14
solver.report = True
solver.max_it = 10000

ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()

#opts[f"{option_prefix}ksp_type"] = "cg" # Extremely slow. Never finishes?
#opts[f"{option_prefix}pc_type"] = "gamg" # Extremely slow. Never finishes?
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "ksp" # Quick execution of the code.
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()

The code runs reasonably quickly with these options, but if I pick “cg” or “gamg”, it takes forever and I don’t know if it’s stuck or making a slow progress (I always end up aborting).