Hello dokken and thanks for helping me that much. On the phone, so hard to make a high quality reply.
I try to copy the relevant information below.
# Define ME function space
deg = 3
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
s_xx, s_yy = seebeck_components
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)
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)
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
if Bridgman_term_bool:
F_T += 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.
Jac = derivative(weak_form,TempVolt,dTV)
# 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-9
solver.report = True
solver.max_it = 10000
ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}pc_type"] = "lu"
opts[f"{option_prefix}pc_type"] = "ksp" # Quick code.
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
opts[f"{option_prefix}ksp_max_it"]= 10000
opts[f"{option_prefix}mat_mumps_icntl_24"] = 1
opts[f"{option_prefix}mat_mumps_icntl_25"] = 0
ksp.setFromOptions()
# Specify the null space, because voltage is defined up to a constant and we don't use any Dirichlet boundary condition to fix the voltage.
nullspace = PETSc.NullSpace().create(constant=True) # type: ignore
PETSc.Mat.setNearNullSpace(solver.A, nullspace)
log.set_log_level(log.LogLevel.INFO)
n, converged = solver.solve(TempVolt)
assert(converged)
print(f'''------- Processing --------
Number of interations: {n:d}''')
# Compute fluxes on boundaries
n = FacetNormal(mesh)
down_heat_flux = form(dot(-κ * grad(temp) + Seebeck_tensor * temp * J_vector, n)*ds(inj_current_curve))
Q_out = assemble_scalar(down_heat_flux)
top_heat_flux = form(dot(-κ * grad(temp) + Seebeck_tensor * temp * J_vector, n)*ds(out_current_curve))
Q_in = assemble_scalar(top_heat_flux)
if thermoelectric_sim:
# Heat flux. (Non TE).
Q = functionspace(mesh, ("CG", deg-1, (mesh.geometry.dim,)))
q = Function(Q)
flux_calculator = Expression(-κ * grad(temp), Q.element.interpolation_points())
q.interpolate(flux_calculator)
# Current density.
J = functionspace(mesh, ("DG", deg-1, (mesh.geometry.dim,)))
Je = Function(J)
Je_flux_calculator = Expression(J_vector, J.element.interpolation_points())
Je.interpolate(Je_flux_calculator)
element = FiniteElement("DG", mesh.ufl_cell(), 1)
J_x = FunctionSpace(mesh, element)
Je_x = Function(J_x)
Je_x_calculator = Expression(J_vector[0].dx(0), J_x.element.interpolation_points())
Je_x.interpolate(Je_x_calculator)
element = FiniteElement("DG", mesh.ufl_cell(), 1)
J_y = FunctionSpace(mesh, element)
Je_y = Function(J_y)
Je_y_calculator = Expression(J_vector[1].dx(1), J_y.element.interpolation_points())
Je_y.interpolate(Je_y_calculator)
# Bridgman heat.
element = FiniteElement("DG", mesh.ufl_cell(), 1)
Q_B = FunctionSpace(mesh, element)
QB = Function(Q_B)
QB_flux_calculator = Expression(- temp * (s_xx * J_vector[0].dx(0) + s_yy * J_vector[1].dx(1)), Q_B.element.interpolation_points(), comm=MPI.COMM_SELF)
QB.interpolate(QB_flux_calculator)
output_0 = TempVolt.sub(0).collapse()
output_1 = TempVolt.sub(1).collapse()
output_2 = q
output_3 = Je
output_4 = QB
output_5 = Je_x
output_6 = Je_y
return output_0, output_1, output_2, output_3, output_4, output_5, output_6
I still have no clue why that term would be any different than the others. In my mind, it does depend on the mesh. It involves a sum of derivatives of a current density which is itself a grad(temp) and grad(volt) quantity, supposedly defined in a DG space of degree 2 if I understand well, since both temp and volt are defined over a FE space of order 3, CG.