Hi,
this works exactly the same, independently of mechanics. If you take the action of the thermal bilinear form on a field with unit value on a given boundary, you will get the total flux on this boundary in a consistent manner.
1 Like
Thanks a lot Jeremy, I am very excited, I will post my (failed) attempt when I have the time. So far my whole personal project solely depends on evaluating these fluxes as accurately as possible, I was thinking about ways to bypass the problem but I did not find any remedy. I now have big hopes!
My current attempt looks like so:
import numpy as np
from dolfinx import log, fem
from dolfinx.fem import (Constant, Function, functionspace, FunctionSpace, assemble_scalar,
dirichletbc, form, locate_dofs_geometrical, VectorFunctionSpace,
locate_dofs_topological, Expression, create_matrix, assemble_matrix)
from dolfinx.fem.petsc import NonlinearProblem
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, curl, as_vector, action)
import gmsh
from dolfinx.io import gmshio, VTXWriter
mesh, cell_markers, facet_markers = gmshio.read_from_msh("meshes/mesh2.msh", MPI.COMM_WORLD, gdim=2)
def thermal_problem(ÎT=0):
current_curves=(1,2)
voltages_curves=(1,2)
inj_current_curve, out_current_curve = current_curves
reading_voltage_curve_0, reading_voltage_curve_1 = voltages_curves
# Define FE function space
deg = 1
el = FiniteElement("CG", mesh.ufl_cell(), deg)
V = FunctionSpace(mesh, el)
u = TestFunction(V)
temp = Function(V)
Îș = 1.8
T_cold = 300
# Define the boundary conditions
left_facets = facet_markers.find(inj_current_curve)
right_facets = facet_markers.find(out_current_curve)
left_dofs_temp = locate_dofs_topological(V, mesh.topology.dim-1, left_facets)
right_dofs_temp = locate_dofs_topological(V, mesh.topology.dim-1, right_facets)
bc_temp_left = dirichletbc(ScalarType(T_cold), left_dofs_temp, V)
bc_temp_right = dirichletbc(ScalarType(T_cold + ÎT), right_dofs_temp, V)
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)
# Weak form.
weak_form = dot(-Îș * grad(temp), grad(u)) * dx
print(f''' ------- Pre-processing --------
Length of the side where heat enters: {assemble_scalar(form(1 * ds(inj_current_curve, domain=mesh)))}
Length of the side where heat leaves the material: {assemble_scalar(form(1 * ds(out_current_curve, domain=mesh)))}
''')
# Solve the PDE.
problem = NonlinearProblem(weak_form, temp, bcs=bcs)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}pc_type"] = "lu"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()
log.set_log_level(log.LogLevel.WARNING)
n, converged = solver.solve(temp)
assert (converged)
print(f'''------- Processing --------
Number of interations: {n:d}''')
avg_temp = assemble_scalar(form(temp * dx(domain=mesh))) / assemble_scalar(form(1 * dx(domain=mesh)))
print(f'avg_temp: {avg_temp}')
# Compute fluxes on boundaries
n = FacetNormal(mesh)
down_heat_flux = form(dot(-Îș * grad(temp), n)*ds(inj_current_curve))
Q_out = assemble_scalar(down_heat_flux)
top_heat_flux = form(dot(-Îș * grad(temp), n)*ds(out_current_curve))
Q_in = assemble_scalar(top_heat_flux)
print(f'''------- Post processing --------
Q_in: {Q_in}
Q_out: {Q_out}''')
Q = functionspace(mesh, ("DG", deg-1, (mesh.geometry.dim,)))
q = Function(Q)
flux_calculator = Expression(-Îș * grad(temp), Q.element.interpolation_points())
q.interpolate(flux_calculator)
# Try an alternate way to compute the fluxes.
virtual_temp = Function(V)
virtual_heat_flux_form = form(action(weak_form, virtual_temp))
virtual_heat_flux = assemble_scalar(virtual_heat_flux_form)
print(f'Virtual heat flux computation: {virtual_heat_flux}')
return Q_in, Q_out
t_b = thermal_problem(ÎT=10)
print(t_b)
Producing the output:
------- Pre-processing --------
Length of the side where heat enters: 0.10000000000000002
Length of the side where heat leaves the material: 0.09999999999999998
------- Processing --------
Number of interations: 1
avg_temp: 305.4445644560497
------- Post processing --------
Q_in: -7.434551004310087
Q_out: 7.494950253817867
Virtual heat flux computation: 0.0
(-7.434551004310087, 7.494950253817867)
In other words, it computes as 0 heat flux instead of something like 7.49 or so.
My mesh is similar to the file.msh I posted in post #10 just above.
Any help is appreciated.
Brief update. Big thank you to you, Jeremy, I could implement your idea in my code. The relevant part is:
# Try an alternate way to compute the fluxes.
u_bc = fem.Function(V)
#bcs = [fem.dirichletbc(u_bc, np.append(left_dofs_temp, left_dofs_temp)) ]
bcs = [fem.dirichletbc(u_bc, right_dofs_temp) ]
virtual_temp = Function(V)
u_bc.interpolate(one)
fem.set_bc(virtual_temp.vector, bcs)
virtual_heat_flux_form = form(action(weak_form, virtual_temp))
virtual_heat_flux = assemble_scalar(virtual_heat_flux_form)
print(f'Virtual heat flux computation: {virtual_heat_flux}')
I do not have any analytical expression but as I increase the order of the basis elements to 3, the values seem to converge to what the above code yields.
Iâm just jumping in at the end, so sorry if this is missing the target, but this is a source that might be useful to you: https://doi.org/10.1115/1.4005187
In my opinion, it is one of those underappreciated gems of papers.
Secondly, you seem to be using the primal form of the heat problem. Given your focus on the fluxes, would a mixed form not be more appropriate? Thatâs probably what Jorgen seems to mean in an earlier response.
3 Likes
Hi Stein and thanks for the information. I wonder if this paper is what inspired Jeremy to write down this part of his nice tutorial.
I am not sure to understand what you mean in your 2nd paragraph. In that MWE I solve a simple heat equation \nabla \cdot (\kappa \nabla T)=0, and I am indeed interested in computing \vec J_Q = -\kappa \nabla T integrated over boundaries.
Do you mean I could have rewritten the weak form as the sum of 2 equations, one involving T and another one involving \vec J_Q? I did not understand this from dokken, I probably missed what he was saying, I suppose.
Interesting⊠actually I will create a new thread about my ârealâ problem, because I already deal with a mixed element space and 2 scalar functions and I am interested in computing some fluxes and I do not understand the results I get.
Yes indeed, itâs a pretty common approach if the fluxes are of primary interest. See Mixed formulation for the Poisson equation â DOLFINx 0.9.0 documentation
1 Like