Hello people,

I’m using a Nonlinear solver for what I think is a linear problem, the reason being that in the future I might want to add non linearities, and this would require many changes in the code if I had chosen a linear solver at first, whereas very little changes with the non linear solver.

I have a code that solves 2 coupled PDEs, with a mixed elements approach. The values seem to make sense, at a highly accurate testings, so I trusted the results, but a recent discussion Solving a system of coupled PDEs in fenicsx - #31 by dokken makes me think my problem is either ill-posed or my problem is non linear. Therefore, I would like to make sure the problem is linear/non linear, as a 1st step to debug.

The coupled PDE’s are:

\nabla \cdot \vec J_e =0 and \nabla \cdot (-\kappa \nabla T) + \rho \vec J_e\cdot \vec J_e-T \left( S_{xx} \frac{\partial J_{e_x}}{\partial x} + S_{yy} \frac{\partial J_{e_y}}{\partial y}\right)=0.

Where \vec J_e = -\sigma \nabla V -\sigma S \nabla T where \sigma =1/\rho is a scalar (float), S is a tensor (2D diagonal matrix in my case, where the 2 non zero elements are floats/scalars).

The boundary conditions are of Dirichlet type for T. I keep 2 “Physical curve” at constant temperature. For V, I employ Neumann type boundary condition, I specify that the voltage between 2 surfaces must be such that a total given electric current must flow.

Fenicsx solves for V and T, the voltage and the temperature. I am essentially solving a thermoelectric heat equation. From these 2 quantities, I can retrieve \vec J_e using the relation above, the heat flux, the electric field, and so on and so forth. I can also compute the total Joule heat in the material, and I confirmed that it matches the expected values to very high accuracy. I tested several cases, and always got meaningful results, to me.

I can show you my full code, but I really wish to stick with what is above. Is the above problem linear?

In any case, here is some code:

```
import numpy as np
from dolfinx import log
from dolfinx.fem import (Constant, Function, FunctionSpace, assemble_scalar,
dirichletbc, form, locate_dofs_geometrical, VectorFunctionSpace,
locate_dofs_topological, Expression)
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.io import XDMFFile, gmshio
from dolfinx.nls.petsc import NewtonSolver
from mpi4py import MPI
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
from ufl import (FiniteElement, Measure, MixedElement, SpatialCoordinate,
TestFunction, TrialFunction, as_tensor, dot, dx, grad, inner,
inv, split, derivative, FacetNormal, div)
from mpi4py import MPI
import gmsh
from dolfinx.io import gmshio, VTXWriter, XDMFFile
mesh, cell_markers, facet_markers = gmshio.read_from_msh("meshes/rectangles.msh", MPI.COMM_WORLD, gdim=2)
# Reminder of the mesh:
#Physical Surface("material", 11) = {2, 1, 3};
#Physical Curve("bottom_left", 12) = {10};
#Physical Curve("top_right", 13) = {6};
inj_current_curve = 12
out_current_curve = 13 # This corresponds to the curve the current leaves the material.
reading_voltage_surface_0 = 12
reading_voltage_surface_1 = 13
# Define ME function space
el = FiniteElement("CG", mesh.ufl_cell(), 4)
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 = 3600e-6
S_yy = S_xx * 10
Seebeck_tensor = as_tensor([[S_xx, 0], [0, S_yy]])
# Define the boundary conditions
left_facets = facet_markers.find(inj_current_curve)
right_facets = facet_markers.find(out_current_curve)
left_dofs = locate_dofs_topological(
ME.sub(1), mesh.topology.dim-1, left_facets)
left_dofs_temp = locate_dofs_topological(
ME.sub(0), mesh.topology.dim-1, left_facets)
right_dofs_temp = locate_dofs_topological(
ME.sub(0), mesh.topology.dim-1, right_facets)
T_cold = 300.0
ΔT = 200
bc_temp_left = dirichletbc(ScalarType(T_cold), left_dofs_temp, ME.sub(0))
bc_temp_right = dirichletbc(ScalarType(T_cold + ΔT), right_dofs_temp, ME.sub(0))
bcs = [bc_temp_left, bc_temp_right]
x = SpatialCoordinate(mesh)
dx = Measure("dx", domain=mesh, subdomain_data=cell_markers)
ds = Measure("ds", domain=mesh, subdomain_data=facet_markers)
resistance = -4.55847847099578
the_current = S_xx * ΔT / (2 * resistance) # Current, in amperes.
the_current = 1e-15#1.0
J = the_current / assemble_scalar(form(1 * ds(inj_current_curve, domain=mesh)))
# 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 PDE.
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
Jac = derivative(weak_form,TempVolt,dTV)
print(f''' ------- Pre-processing --------
Current density: {J}
Length of the side where current is injected: {assemble_scalar(form(1 * ds(inj_current_curve, domain=mesh)))}
Length of the side where current leaves the wire: {assemble_scalar(form(1 * ds(out_current_curve, domain=mesh)))}
This should correspond to the current injected: {assemble_scalar(form(J*ds(out_current_curve)))}
''')
# Solve the PDE.
problem = NonlinearProblem(weak_form, TempVolt, bcs=bcs, J = Jac)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-10
solver.report = True
solver.max_it = 10000
ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "ksp"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()
log.set_log_level(log.LogLevel.WARNING)
n, converged = solver.solve(TempVolt)
assert(converged)
print(f'''------- Processing --------
Number of interations: {n:d}''')
avg_voltage_curve_integral_0 = assemble_scalar(form(volt * ds(reading_voltage_surface_0, domain=mesh))) / assemble_scalar(form(1 * ds(reading_voltage_surface_0, domain=mesh)))
avg_voltage_curve_integral_1 = assemble_scalar(form(volt * ds(reading_voltage_surface_1, domain=mesh))) / assemble_scalar(form(1 * ds(reading_voltage_surface_1, domain=mesh)))
voltage_diff = avg_voltage_curve_integral_1 - avg_voltage_curve_integral_0
avg_temperature_surface_integral_2 = assemble_scalar(form(temp * ds(inj_current_curve, domain=mesh))) / assemble_scalar(form(1 * ds(inj_current_curve, domain=mesh)))
avg_temperature_surface_integral_3 = assemble_scalar(form(temp * ds(out_current_curve, domain=mesh))) / assemble_scalar(form(1 * ds(out_current_curve, domain=mesh)))
# Compute fluxes on boundaries
n = FacetNormal(mesh)
down_heat_flux = form(dot(-κ * grad(temp) + Seebeck_tensor * temp * J_vector, n)*ds(out_current_curve))
Q_out = assemble_scalar(down_heat_flux)
top_heat_flux = form(dot(-κ * grad(temp) + Seebeck_tensor * temp * J_vector, n)*ds(inj_current_curve))
Q_in = assemble_scalar(top_heat_flux)
print(f'''--------- Post processing -------
voltage measured: {voltage_diff} V.
voltage left: {avg_voltage_curve_integral_0} V.
voltage right: {avg_voltage_curve_integral_1} V.
the_current: {the_current} A.
This should correspond to the current injected: {assemble_scalar(form(dot(J_vector, n)*ds(out_current_curve)))}
Temperature of current injection contact: {avg_temperature_surface_integral_2} K.
Temperature of current outjection contact: {avg_temperature_surface_integral_3} K.
Computed generated Joule heat through the whole volume: {assemble_scalar(form(dot(rho * J_vector, J_vector) * dx))} W.
Computed Joule heat with I^2R: {the_current**2 * resistance} W.
Make sure div(J_e) = 0 throughought the whole volume: {assemble_scalar(form((J_vector[0].dx(0) + J_vector[1].dx(1)) * dx))}
Make sure div(J_e) = 0 throughought the whole volume: {assemble_scalar(form(div(J_vector) * dx))}
Average temperature of whole volume: {assemble_scalar(form(temp * dx(domain=mesh))) / assemble_scalar(form(1 * dx(domain=mesh)))}
Electrical resistance should be slightly lower than {0.5 /0.1 * rho}. It is worth {voltage_diff / the_current}.
More accurately, it is worth {voltage_diff / the_current + S_xx * ΔT / the_current}
Q_out: {Q_out} W/m^2.
Q_in: {Q_in} W/m^2.
Power produced by the TE element: {Q_in + Q_out} W.
# Heat flux.
Q = VectorFunctionSpace(mesh, ("DG", 3))
q = Function(Q)
flux_calculator = Expression(-κ * grad(temp), Q.element.interpolation_points())
q.interpolate(flux_calculator)
# Current density.
J = VectorFunctionSpace(mesh, ("DG", 3))
Je = Function(J)
Je_flux_calculator = Expression(-sigma * grad(volt) - sigma * Seebeck_tensor * grad(temp), J.element.interpolation_points())
Je.interpolate(Je_flux_calculator)
Output_0 = TempVolt.sub(0).collapse()
Output_1 = TempVolt.sub(1).collapse()
output_3 = q
output_4 = Je
with VTXWriter(MPI.COMM_WORLD, "results/temperature.bp", [Output_0], engine="BP4") as vtx:
vtx.write(0.0)
with VTXWriter(MPI.COMM_WORLD, "results/voltage.bp", [Output_1], engine="BP4") as vtx:
vtx.write(0.0)
with VTXWriter(MPI.COMM_WORLD, "results/heat_flux.bp", [output_3], engine="BP4") as vtx:
vtx.write(0.0)
with VTXWriter(MPI.COMM_WORLD, "results/current_density.bp", [output_4], engine="BP4") as vtx:
vtx.write(0.0)
```