# Computing fluxes accurately

I have a few questions about understanding FEM.
I have noticed that if I define a function space like so: `V = FunctionSpace(mesh, ("CG", 1))`, then I define a variable I want to solve for, in this case the temperature for the heat equation: `Temp = Function(V)`, from which I compute the flux (in this case the heat flux) like so:

``````heat_flux = VectorFunctionSpace(mesh, ("DG", 1))
JQ = Function(heat_flux)
``````

then there will be a “big” inaccuracy in the calculated flux. In my particular case, I have a rod in which I keep a surface at constant temperature (Dirichlet boundary condition), and it loses heat via convection and radiation. In the steady state, i.e. after a large time, the result converges towards the steady state solution, in which the temperature does not change much with respect to time. When I integrate the flux over the surfaces where I keep the temperature fixed, and over all the other surfaces, I should get the same magnitude (with a sign flip). Physically this means that all heat that enters (in W), equals the one that leaves the rod, and the temperature doesn’t evolve anymore.
However I do not get this with the above code. I get that the temperature indeed converges to a good value, but that the outward flux is only about 73% of the inward flux. I didn’t know why, but I tried to modify the above line to `V = FunctionSpace(mesh, ("CG", 2))`. The computations take a significant bigger amount of time to run, but this achieves the correct physical result of inward flux = outward flux in steady state. This is very unfortunate because the simulations takes way longer now, while the temperature was correctly (I believe?) computed in the first place. Only the flux was not.

Is this expected behavior? I.e. whenever we are interested in the flux, we should use at least a 2nd order function space? Is my “fix” a correct way to deal with this issue? Or are there more efficient way to achieve this, e.g. with Paraview?

Without supplying a minimal working example of how you have implemented this, its hard to Give you guidance as to what goes wrong in your code.

I didn’t think it was necessary, but here it goes. The mesh is generated with gmsh, here’s the file.geo:

``````SetFactory("OpenCASCADE");
Box(1) = {0, 0, 0, 1e-3, 1e-3, 1e-2};
Physical Volume("the_rod", 20) = {1};
Physical Surface("left_side",21) = {6};
Physical Surface("res_of_rod_surface",22) = {1,2,3,4,5};
``````

and the corresponding file.msh:

``````\$MeshFormat
4.1 0 8
\$EndMeshFormat
\$PhysicalNames
3
2 21 "left_side"
2 22 "res_of_rod_surface"
3 20 "the_rod"
\$EndPhysicalNames
\$Entities
8 12 6 1
1 0 0 0.01 0
2 0 0 0 0
3 0 0.001 0.01 0
4 0 0.001 0 0
5 0.001 0 0.01 0
6 0.001 0 0 0
7 0.001 0.001 0.01 0
8 0.001 0.001 0 0
1 -1e-07 -1e-07 -1.000000000002735e-07 1e-07 1e-07 0.0100001 0 2 2 -1
2 -1e-07 -9.999999999994822e-08 0.009999900000000001 1e-07 0.0010001 0.0100001 0 2 1 -3
3 -1e-07 0.0009999 -1.000000000002735e-07 1e-07 0.0010001 0.0100001 0 2 4 -3
4 -1e-07 -9.999999999994822e-08 -1e-07 1e-07 0.0010001 1e-07 0 2 2 -4
5 0.0009999 -1e-07 -1.000000000002735e-07 0.0010001 1e-07 0.0100001 0 2 6 -5
6 0.0009999 -9.999999999994822e-08 0.009999900000000001 0.0010001 0.0010001 0.0100001 0 2 5 -7
7 0.0009999 0.0009999 -1.000000000002735e-07 0.0010001 0.0010001 0.0100001 0 2 8 -7
8 0.0009999 -9.999999999994822e-08 -1e-07 0.0010001 0.0010001 1e-07 0 2 6 -8
9 -9.999999999994822e-08 -1e-07 -1e-07 0.0010001 1e-07 1e-07 0 2 2 -6
10 -9.999999999994822e-08 -1e-07 0.009999900000000001 0.0010001 1e-07 0.0100001 0 2 1 -5
11 -9.999999999994822e-08 0.0009999 -1e-07 0.0010001 0.0010001 1e-07 0 2 4 -8
12 -9.999999999994822e-08 0.0009999 0.009999900000000001 0.0010001 0.0010001 0.0100001 0 2 3 -7
1 -1e-07 -9.999999999994822e-08 -1.000000000002735e-07 1e-07 0.0010001 0.0100001 1 22 4 -1 4 3 -2
2 0.0009999 -9.999999999994822e-08 -1.000000000002735e-07 0.0010001 0.0010001 0.0100001 1 22 4 -5 8 7 -6
3 -9.999999999994822e-08 -1e-07 -1.000000000002735e-07 0.0010001 1e-07 0.0100001 1 22 4 -9 1 10 -5
4 -9.999999999994822e-08 0.0009999 -1.000000000002735e-07 0.0010001 0.0010001 0.0100001 1 22 4 -11 3 12 -7
5 -9.999999999994822e-08 -9.999999999994822e-08 -1e-07 0.0010001 0.0010001 1e-07 1 22 4 -4 9 8 -11
6 -9.999999999994822e-08 -9.999999999994822e-08 0.009999900000000001 0.0010001 0.0010001 0.0100001 1 21 4 -2 10 6 -12
1 -9.999999999994822e-08 -9.999999999994822e-08 -1.000000000002735e-07 0.0010001 0.0010001 0.0100001 1 20 6 -1 2 -3 4 -5 6
\$EndEntities
\$Nodes
19 86 1 86
0 1 0 1
1
0 0 0.01
0 2 0 1
2
0 0 0
0 3 0 1
3
0 0.001 0.01
0 4 0 1
4
0 0.001 0
0 5 0 1
5
0.001 0 0.01
0 6 0 1
6
0.001 0 0
0 7 0 1
7
0.001 0.001 0.01
0 8 0 1
8
0.001 0.001 0
1 1 0 9
9
10
11
12
13
14
15
16
17
0 0 0.001000000000000003
0 0 0.002000000000000007
0 0 0.00300000000000001
0 0 0.004000000000000013
0 0 0.005000000000000011
0 0 0.006000000000000008
0 0 0.007000000000000006
0 0 0.008000000000000004
0 0 0.009000000000000003
1 3 0 9
18
19
20
21
22
23
24
25
26
0 0.001 0.001000000000000003
0 0.001 0.002000000000000007
0 0.001 0.00300000000000001
0 0.001 0.004000000000000013
0 0.001 0.005000000000000011
0 0.001 0.006000000000000008
0 0.001 0.007000000000000006
0 0.001 0.008000000000000004
0 0.001 0.009000000000000003
1 5 0 9
27
28
29
30
31
32
33
34
35
0.001 0 0.001000000000000003
0.001 0 0.002000000000000007
0.001 0 0.00300000000000001
0.001 0 0.004000000000000013
0.001 0 0.005000000000000011
0.001 0 0.006000000000000008
0.001 0 0.007000000000000006
0.001 0 0.008000000000000004
0.001 0 0.009000000000000003
1 7 0 9
36
37
38
39
40
41
42
43
44
0.001 0.001 0.001000000000000003
0.001 0.001 0.002000000000000007
0.001 0.001 0.00300000000000001
0.001 0.001 0.004000000000000013
0.001 0.001 0.005000000000000011
0.001 0.001 0.006000000000000008
0.001 0.001 0.007000000000000006
0.001 0.001 0.008000000000000004
0.001 0.001 0.009000000000000003
2 1 0 10
45
46
47
48
49
50
51
52
53
54
0 0.0005 0.003500000000000011
0 0.0005 0.006500000000000007
0 0.0005 0.009500000000000001
0 0.0005 0.0005000000000000016
0 0.0005 0.001500000000000005
0 0.0005 0.005500000000000008
0 0.0004999999999999999 0.002500000000000008
0 0.0005000000000000008 0.004500000000000013
0 0.0005000000000000008 0.007500000000000008
0 0.0005 0.008499999999999997
2 2 0 10
55
56
57
58
59
60
61
62
63
64
0.001 0.0005 0.003500000000000011
0.001 0.0005 0.006500000000000007
0.001 0.0005 0.009500000000000001
0.001 0.0005 0.0005000000000000016
0.001 0.0004999999999999999 0.002500000000000008
0.001 0.0005 0.004500000000000013
0.001 0.0005 0.001500000000000005
0.001 0.0005 0.005500000000000009
0.001 0.0005000000000000008 0.007500000000000008
0.001 0.0005 0.008499999999999997
2 3 0 10
65
66
67
68
69
70
71
72
73
74
0.0005000000000000016 0 0.008499999999999997
0.0005 0 0.003500000000000011
0.0005 0 0.006500000000000007
0.0005 0 0.0005000000000000016
0.0004999999999999999 0 0.002500000000000008
0.0005000000000000008 0 0.004500000000000013
0.0005 0 0.001500000000000005
0.0005 0 0.005500000000000008
0.0005 0 0.007500000000000008
0.0005 0 0.009500000000000001
2 4 0 10
75
76
77
78
79
80
81
82
83
84
0.0005000000000000016 0.001 0.008499999999999997
0.0005 0.001 0.003500000000000011
0.0004999999999999999 0.001 0.006500000000000008
0.0004999999999999999 0.001 0.0005000000000000016
0.0005 0.001 0.001500000000000005
0.0004999999999999999 0.001 0.002500000000000008
0.0005000000000000008 0.001 0.004500000000000013
0.0005000000000000001 0.001 0.005500000000000008
0.0005000000000000003 0.001 0.007500000000000005
0.0005 0.001 0.009500000000000001
2 5 0 1
85
0.0005 0.0005 0
2 6 0 1
86
0.0005 0.0005 0.01
3 1 0 0
\$EndNodes
\$Elements
7 366 1 366
2 1 2 40
1 1 3 47
2 17 1 47
3 4 2 48
4 2 9 48
5 3 26 47
6 18 4 48
7 9 10 49
8 9 18 48
9 18 9 49
10 10 11 51
11 10 19 49
12 19 10 51
13 11 12 45
14 20 11 45
15 11 20 51
16 12 13 52
17 12 21 45
18 21 12 52
19 13 14 50
20 13 50 52
21 14 15 46
22 14 46 50
23 15 16 53
24 15 24 46
25 24 15 53
26 16 17 54
27 53 16 54
28 26 17 47
29 17 26 54
30 19 18 49
31 20 19 51
32 21 20 45
33 22 21 52
34 23 22 50
35 50 22 52
36 24 23 46
37 46 23 50
38 25 24 53
39 26 25 54
40 25 53 54
2 2 2 40
41 5 57 7
42 35 57 5
43 8 58 6
44 6 58 27
45 7 57 44
46 36 58 8
47 27 61 28
48 27 58 36
49 36 61 27
50 28 59 29
51 37 59 28
52 28 61 37
53 29 55 30
54 38 55 29
55 29 59 38
56 30 60 31
57 30 55 39
58 39 60 30
59 31 62 32
60 60 62 31
61 32 56 33
62 32 62 56
63 33 63 34
64 33 56 42
65 42 63 33
66 34 64 35
67 63 64 34
68 44 57 35
69 35 64 44
70 37 61 36
71 38 59 37
72 39 55 38
73 40 60 39
74 41 62 40
75 40 62 60
76 42 56 41
77 56 62 41
78 43 63 42
79 44 64 43
80 43 64 63
2 3 2 40
81 5 1 74
82 1 17 74
83 2 6 68
84 9 2 68
85 35 5 74
86 6 27 68
87 10 9 71
88 27 9 68
89 9 27 71
90 11 10 69
91 10 28 69
92 28 10 71
93 12 11 66
94 11 29 66
95 29 11 69
96 13 12 70
97 30 12 66
98 12 30 70
99 14 13 72
100 13 70 72
101 15 14 67
102 67 14 72
103 16 15 73
104 33 15 67
105 15 33 73
106 17 16 65
107 65 16 73
108 35 17 65
109 17 35 74
110 27 28 71
111 28 29 69
112 29 30 66
113 30 31 70
114 31 32 72
115 70 31 72
116 32 33 67
117 32 67 72
118 33 34 73
119 34 35 65
120 34 65 73
2 4 2 40
121 7 84 3
122 3 84 26
123 4 78 8
124 18 78 4
125 44 84 7
126 8 78 36
127 19 79 18
128 36 78 18
129 18 79 36
130 20 80 19
131 37 79 19
132 19 80 37
133 21 76 20
134 20 76 38
135 38 80 20
136 22 81 21
137 39 76 21
138 21 81 39
139 23 82 22
140 22 82 81
141 24 77 23
142 77 82 23
143 25 83 24
144 24 83 77
145 26 75 25
146 75 83 25
147 44 75 26
148 26 84 44
149 36 79 37
150 37 80 38
151 38 76 39
152 39 81 40
153 40 82 41
154 81 82 40
155 41 77 42
156 41 82 77
157 42 83 43
158 77 83 42
159 43 75 44
160 43 83 75
2 5 2 4
161 2 4 85
162 6 2 85
163 4 8 85
164 8 6 85
2 6 2 4
165 1 86 3
166 5 86 1
167 3 86 7
168 7 86 5
3 1 4 198
169 74 47 84 44
170 84 57 74 44
171 68 61 36 79
172 79 9 48 68
173 9 61 68 79
174 36 68 79 48
175 69 51 59 28
176 80 51 19 59
177 36 48 78 58
178 38 76 45 55
179 36 68 48 58
180 55 45 11 66
181 59 38 51 69
182 10 49 79 71
183 79 49 9 71
184 61 10 79 71
185 9 61 79 71
186 51 80 38 59
187 39 45 76 55
188 12 55 66 45
189 47 26 84 44
190 74 17 47 35
191 57 35 74 44
192 47 74 35 44
193 27 9 61 68
194 36 18 48 79
195 27 36 68 61
196 9 18 79 48
197 36 27 68 58
198 69 10 51 28
199 37 80 19 59
200 79 49 10 19
201 76 45 20 38
202 61 10 71 28
203 59 38 69 29
204 11 55 66 29
205 18 36 48 78
206 18 9 79 49
207 11 51 38 69
208 10 79 19 61
209 37 61 19 79
210 45 11 38 55
211 9 27 61 71
212 12 45 39 55
213 59 51 19 28
214 30 12 55 66
215 21 45 76 39
216 80 20 51 38
217 39 30 70 55
218 12 30 55 70
219 26 17 35 47
220 26 35 44 47
221 39 12 55 70
222 51 10 19 28
223 61 19 10 28
224 28 37 61 19
225 69 38 11 29
226 28 37 19 59
227 55 11 38 29
228 11 20 38 51
229 20 11 38 45
230 12 21 39 45
231 54 53 83 16
232 73 53 16 83
233 64 54 43 65
234 35 65 54 64
235 54 64 43 75
236 54 64 75 44
237 65 73 16 63
238 64 63 65 43
239 83 16 73 63
240 54 17 65 35
241 54 35 64 44
242 43 54 75 83
243 75 26 54 44
244 83 77 42 63
245 81 60 22 70
246 33 42 77 63
247 81 22 52 70
248 52 72 70 22
249 83 73 33 63
250 15 53 73 83
251 62 56 41 72
252 65 16 54 83
253 26 35 54 44
254 17 26 35 54
255 22 82 72 62
256 52 50 72 22
257 82 81 60 22
258 82 60 62 22
259 60 81 39 70
260 43 54 83 65
261 39 81 52 70
262 82 50 22 72
263 72 82 41 62
264 14 46 67 56
265 82 72 41 50
266 83 65 16 63
267 86 57 74 84
268 22 72 70 60
269 15 33 83 73
270 83 43 65 63
271 72 67 56 14
272 82 46 50 41
273 56 50 72 14
274 22 72 60 62
275 46 56 33 67
276 74 47 86 84
277 50 56 72 41
278 56 46 33 77
279 41 77 46 56
280 33 77 46 83
281 50 56 46 14
282 50 56 41 46
283 53 15 24 83
284 46 33 83 15
285 62 56 72 32
286 46 33 15 67
287 82 77 46 41
288 39 12 70 52
289 30 39 70 60
290 42 33 77 56
291 33 77 83 63
292 52 21 81 39
293 72 67 32 56
294 46 77 24 83
295 77 82 46 23
296 54 53 25 83
297 24 15 46 83
298 82 46 23 50
299 52 72 13 70
300 58 48 85 68
301 54 25 75 83
302 78 48 85 58
303 60 31 72 70
304 21 12 39 52
305 81 82 60 40
306 50 52 72 13
307 73 65 34 63
308 62 31 72 60
309 82 60 40 62
310 64 63 34 65
311 74 1 5 86
312 47 3 1 86
313 19 79 49 18
314 78 4 48 18
315 68 6 58 27
316 37 79 36 61
317 9 2 48 68
318 49 10 9 71
319 45 12 11 66
320 37 38 80 59
321 36 78 8 58
322 21 45 20 76
323 76 38 39 55
324 19 20 51 80
325 30 66 55 29
326 51 11 10 69
327 27 61 71 28
328 5 7 57 86
329 86 3 7 84
330 29 69 59 28
331 85 6 8 58
332 85 4 2 48
333 17 47 1 74
334 35 74 5 57
335 47 26 3 84
336 84 7 57 44
337 85 2 6 68
338 85 8 4 78
339 57 74 5 86
340 74 47 1 86
341 78 85 8 58
342 2 85 48 68
343 58 85 6 68
344 48 85 4 78
345 86 47 3 84
346 7 57 86 84
347 34 65 35 64
348 43 75 64 44
349 54 25 26 75
350 17 16 54 65
351 16 15 53 73
352 34 33 73 63
353 83 42 43 63
354 24 25 53 83
355 41 42 77 56
356 46 15 14 67
357 32 31 72 62
358 60 31 70 30
359 77 24 23 46
360 14 13 50 72
361 40 41 82 62
362 22 52 21 81
363 22 50 82 23
364 40 60 81 39
365 33 32 67 56
366 52 13 12 70
\$EndElements
``````

Now the FEniCSx code:

``````from dolfinx.fem import (Constant, dirichletbc, Function, FunctionSpace, assemble_scalar, Expression,
form, locate_dofs_geometrical, locate_dofs_topological, VectorFunctionSpace)
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.io import gmshio, XDMFFile
from dolfinx.nls.petsc import NewtonSolver
import gmsh
import meshio
from mpi4py import MPI
import numpy as np
from petsc4py.PETSc import ScalarType, Options
from ufl import (SpatialCoordinate, TestFunction, Measure, dot,
dx, grad, inner, MixedElement, FiniteElement, FacetNormal)
import matplotlib.pyplot as plt

proc = MPI.COMM_WORLD.rank

def create_mesh(mesh, cell_type, prune_z=False):
cells = mesh.get_cells_type(cell_type)
cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
points = mesh.points[:,:2] if prune_z else mesh.points
out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
return out_mesh

mesh, cell_markers, facet_markers = gmshio.read_from_msh("the_mesh.msh", MPI.COMM_WORLD, gdim=3)
if proc == 0:
triangle_mesh = create_mesh(msh, "triangle", prune_z=False)
tetra_mesh = create_mesh(msh, "tetra", prune_z=False)
meshio.write('tetra.xdmf', tetra_mesh)
meshio.write('mt.xdmf', triangle_mesh)

with XDMFFile(MPI.COMM_WORLD, 'tetra.xdmf', "r") as xdmf:

mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim-1)
with XDMFFile(MPI.COMM_WORLD, 'mt.xdmf', "r") as xdmf:

# Define function space
V = FunctionSpace(mesh, ("CG", 1))

# Define the "unknown" variable of the equation.
Temp = Function(V)

heat_flux = VectorFunctionSpace(mesh, ("DG", 1))
JQ = Function(heat_flux)

# Physical parameters.
κ = 49.0
C_p = 120
density = 9.7e3 # kg / m^3.
T_amb = 0.0
σ_SB = 5.6704e-8
h = 10

T_fixed_surface = 21
convection_surface = 22
T_left = 100

# Define Dirichlet boundary conditions to the facets corresponding to the left end of the rod.
left_facets = ft.find(T_fixed_surface)
left_dofs = locate_dofs_topological(V, mesh.topology.dim-1, left_facets)
bc = dirichletbc(ScalarType(T_left), left_dofs, V)
bcs = [bc]

x = SpatialCoordinate(mesh)
dx = Measure("dx", domain=mesh,subdomain_data=ct)
ds = Measure("ds", domain=mesh, subdomain_data=ft)

# Initial temperature of the rod.
T_0 = 0.0

def temp_init(x):
values = np.full(x.shape[1], T_0, dtype = ScalarType)
return values

Temp.name = "Temp"
Temp.interpolate(temp_init)

# Simulation parameters.
t = 0 # Start time
t_final = 2.0 # Final time
num_steps = 260
dt = t_final / num_steps # time step size

T_n = Function(V)
T_n.name = "T_n"
T_n.interpolate(temp_init)

v = TestFunction(V)

F_T = density * C_p * (Temp - T_n )* v * dx + dt * dot(κ * grad(Temp), grad(v)) * dx + v * h * (Temp - T_amb) * ds(convection_surface)

# Define the problem.
problem = NonlinearProblem(F_T, Temp, bcs=bcs)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-14
solver.report = True

ksp = solver.krylov_solver
opts = Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "cg"
opts[f"{option_prefix}pc_type"] = "gamg"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
opts[f"{option_prefix}ksp_max_it"]= 10000

ksp.setFromOptions()

times, Qs_ratios, spatial_average_temperatures = [], [], []
for time_step in np.linspace(t, t_final, num_steps):
times.append(t)
# Compute the heat flux.
JQ_expr = Expression(-κ * grad(Temp), heat_flux.element.interpolation_points())
JQ.interpolate(JQ_expr)

# Compute total heat convected away.
n = FacetNormal(mesh)

# Heat entering the rod from the left side.
Q_in = assemble_scalar(form(dot(JQ, n)*ds(T_fixed_surface)))

# Heat leaving the rod on all other surfaces.
Q_out = assemble_scalar(form(dot(JQ, n)*ds(convection_surface)))

if Q_in == 0:
Qs_ratio = 0
else:
Qs_ratio = abs(Q_out) / abs(Q_in)
Qs_ratios.append(Qs_ratio * 100)
t += dt

n, converged = solver.solve(Temp)

T_n.x.array[:] = Temp.x.array

spatial_average_temperature = assemble_scalar(form(Temp * dx(domain = mesh))) / assemble_scalar(form(1 * dx(domain=mesh)))
spatial_average_temperatures.append(spatial_average_temperature)

plt.plot(times, spatial_average_temperatures)
plt.plot(times, Qs_ratios)
plt.show()
plt.close()
``````

This produces a solution to the heat equation that reaches steady state in about 0.5 s. However the ratio of the flux entering and the flux leaving the rod is about 0.195, which doesn’t make physical sense. If this was really true, then the rod would be heating at a high rate, whereas the correct computation of the temperature profile shows that temperature has stabilized quickly. The average temperature is near 30.5.
But if I do the single line modification, i.e. increase the degree of the function space to 2, the flux ratio reaches 0.945 which is much better (in my real code, it is closer to 100). Average T is near 30.5, same as before.

Graphically, you can see the behavior. Blue curve is average temperature, orange curve is the fluxes ratio which should converge to 100 (it’s in %).

vs

First, it seems like you have a scaling mistake in your code:

What happened to `dt` for the `ds` term.

Secondly, your code can be simplified quite alot:

``````mesh, cell_markers, facet_markers = gmshio.read_from_msh("the_mesh.msh", MPI.COMM_WORLD, gdim=3)
# Define function space
V = FunctionSpace(mesh, ("Lagrange", 2))
# Define the "unknown" variable of the equation.
Temp = Function(V)

# Physical parameters.
κ = 49.0
C_p = 120
density = 9.7e3 # kg / m^3.
T_amb = 0.0
q_SB = 5.6704e-8
h = 10

T_fixed_surface = 21
convection_surface = 22
T_left = 100

# Define Dirichlet boundary conditions to the facets corresponding to the left end of the rod.
left_facets = facet_markers.find(T_fixed_surface)
left_dofs = locate_dofs_topological(V, mesh.topology.dim-1, left_facets)
bc = dirichletbc(ScalarType(T_left), left_dofs, V)
bcs = [bc]

x = SpatialCoordinate(mesh)
dx = Measure("dx", domain=mesh,subdomain_data=cell_markers)
ds = Measure("ds", domain=mesh, subdomain_data=facet_markers)

# Initial temperature of the rod.
T_0 = 0.0

def temp_init(x):
values = np.full(x.shape[1], T_0, dtype = ScalarType)
return values

Temp.name = "Temp"
Temp.interpolate(temp_init)

# Simulation parameters.
t = 0 # Start time
t_final = 6.0 # Final time
num_steps = 100
dt = t_final / num_steps # time step size

T_n = Function(V)
T_n.name = "T_n"
T_n.interpolate(temp_init)

v = TestFunction(V)

F_T = density * C_p * (Temp - T_n )* v * dx + dt * dot(κ * grad(Temp), grad(v)) * dx + dt * v * h * (Temp - T_amb) * ds(convection_surface)

# Define the problem.
problem = NonlinearProblem(F_T, Temp, bcs=bcs)
solver = NewtonSolver(mesh.comm, problem)
solver.convergence_criterion = "residual"
solver.rtol = 1e-14
solver.report = True

ksp = solver.krylov_solver
opts = Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
opts[f"{option_prefix}ksp_max_it"]= 10000

ksp.setFromOptions()

times, Qs_ratios, spatial_average_temperatures = [], [], []
n = FacetNormal(mesh)
left_flux = form(dot(-κ * grad(Temp), n)*ds(T_fixed_surface))
right_flux = form(dot(-κ * grad(Temp), n)*ds(convection_surface))

Q = VectorFunctionSpace(mesh, ("DG", 1))
q = Function(Q)
flux_calculator = Expression(-κ * grad(Temp), Q.element.interpolation_points())

out_file = VTXWriter(mesh.comm, "t.bp", [Temp], engine="BP4")
out_flux = VTXWriter(mesh.comm, "flux.bp", [q], engine="BP4")

#set_log_level(LogLevel.INFO)
for time_step in np.linspace(t, t_final, num_steps):
times.append(t)

# Compute total heat convected away.

# Heat entering the rod from the left side.
Q_in = assemble_scalar(left_flux)

# Heat leaving the rod on all other surfaces.
Q_out = assemble_scalar(right_flux)

if Q_in == 0:
Qs_ratio = 0
else:
Qs_ratio = abs(Q_out) / abs(Q_in)
print(Q_out, Q_in, time_step)
Qs_ratios.append(Qs_ratio * 100)
t += dt
out_file.write(t)
n, converged = solver.solve(Temp)

T_n.x.array[:] = Temp.x.array
q.interpolate(flux_calculator)
out_flux.write(time_step)
#spatial_average_temperature = assemble_scalar(form(Temp * dx(domain = mesh))) / assemble_scalar(form(1 * dx(domain=mesh)))
#spatial_average_temperatures.append(spatial_average_temperature)

out_flux.close()
out_file.close()
#plt.plot(times, spatial_average_temperatures)
plt.plot(times, Qs_ratios)
plt.savefig("Test.png")
plt.close()
``````

I would also like to note that your mesh is incredibly coarse (and it is scaled such that the cell jacobian is very tiny, which could lead to numerical issues. For instance the fluxes are very large (1e6 magnitude)

You are also not attaching the appropriate null space (the constant space) to your multi-grid solver.

Considering the code above with P1, P2, P3 gives

First of all, I really thank you dokken, I did not expect so much valuable feedback on my code. I really appreciate the simplification regarding the mesh. Thanks for pointing out the scaling issue with the dt, it is also missing in my code (although not my old code, thankfully!).
I chose a coarse mesh for the MWE, because I have to paste the full file.msh here, and there is a characters limit in this discourse website, so I cannot paste my real mesh(es), this happened in the past also.

Could you please elaborate what you mean by

You are also not attaching the appropriate null space (the constant space) to your multi-grid solver.

I honestly have no idea what this means, and I don’t know if this means my results are wrong due to this, and I also don’t know if you fixed it with the code you presented here.

And by

Considering the code above with P1, P2, P3 gives

Do you mean you tried the code and changed the line `V = FunctionSpace(mesh, ("Lagrange", 2))`, modifying the 2 by a 1 and by a 3? Or something else?

I have ran your code (with a small modification, “engine” was not recognized, I guess I have an older dolfinx version, which was installed with anaconda on this machine). The behavior I describe in my original post holds. The flux ratio is 99.1796898933635 when I chose a final time of 20 s, 450 time steps, and , `V = FunctionSpace(mesh, ("CG", 3))`, `Q = VectorFunctionSpace(mesh, ("DG", 2))`. This makes physical sense. However if I specify degrees of 1 for V and Q, the value is only worth 6.816261089903455! This doesn’t make any physical sense…

So I still have my original question. If I look at your 3 pictures, I see no difference between the cases, I do not really understand how to reconciliates with the behavior I observe regarding Qratios.

When using iterative solvers such as GAMG, it is often important to consider the near null space of the operator `A`. However, with your current implementation of a diffusion equation, I don’t believe it is an issue of great concern.

Yes, changing the degree of the space.

I would say there is a quite massive difference (especially when choosing a coarse mesh).

The fluxes in the first order case are constant per element (as you can clearly see by the first plot above), while in the case of second or third order Lagrange, the fluxes are in DG-1/DG-2 for 2/3. Thus they vary across elements, yielding a huge difference where there is rapid change in the diffusion (which then is integrated and accumulated).

I see, thanks for all the information.
I didn’t see the difference of colors in the plots, but now that you mention it, I see it. Still surprisingly perhaps, the value of the range doesn’t seem impacted (Paraview’s scale is the same for all 3 plots). But the computations over the surfaces are.

Therefore, can I conclude that if I want to compute fluxes accurately, I need a higher than 1 degree Function space for the variable from which the flux is computed? Or is there another way to deduce the flux more accurately without consuming as much CPU power?

Hello again @dokken , I have a feeling my code computes fluxes innacurately on some meshes (even 2D meshes). I believe the reason is that for corners, n, the facetnormal, might be badly evaluated or not even defined. As a result, the flux near or exactly at these points is enormous in magnitude.

My big problem is that unfortunately this isn’t just a cosmetic visual problem if the heat flux is badly computed, because I am integrating the flux over a surface in order to compute the total heat and use it elsewhere, so I absolutely need to compute the flux accurately. Are there tricks to avoid these problematic evaluations?

If you want to avoid issues with inacurate normals, have a look at: Setting scalar function for BDM boundary condition with controlled flux [mixed Poisson] - #9 by dokken which can evaluate the FacetNormals exactly for every cell.

As you have not provided a new minimal example to illustrate how you currently want to use the fluxes with another program, I cannot give any further guidance.

1 Like

Thanks a lot @dokken. Actually I am not 100% sure this can fix my “problem”, but I think the facetnormals have an issue with corners in my meshes. In some cases, I have an “L” shaped mesh, where heat enters and leave at the ends, and I wish to have accurate fluxes on those curves. Near the right angle in the “L” the flux is also problematic, but I think this matters less, as I do not compute it over a curve.

I have checked your code, and it’s somewhat above my head, I also noticed it has been included in a “test” file of dolfinx? So, not usable as a classical import? I haven’t been able to use it with my mesh.

Here is the mesh file:

``````\$MeshFormat
4.1 0 8
\$EndMeshFormat
\$PhysicalNames
3
1 1 "bottom_left"
1 2 "top_right"
2 3 "material"
\$EndPhysicalNames
\$Entities
12 4 1 0
1 0 0 0 0
2 0.1 0 0 0
3 0.3 0.2 0 0
4 0.3 0.3 0 0
5 0.1716932264052801 0.1095292311834418 0 0
6 0.2303608967257341 0.590084103381926 0 0
7 0.2889801197747777 0.469819750698045 0 0
8 0.2942944418450887 0.4724398881018764 0 0
9 0.1970192363160563 0.06675102825907428 0 0
10 0.2154768708222802 0.07208204727133675 0 0
11 0.2693326138716083 0.1789394416980217 0 0
12 0.2732394432029431 0.2925042062925737 0 0
1 -1.000000000028756e-07 -1e-07 -1e-07 0.1000001 1e-07 1e-07 1 1 2 1 -2
2 0.2999999 0.1999999 -1e-07 0.3000001 0.3000001 1e-07 1 2 2 3 -4
3 -1.000000000028756e-07 -1.000000000028756e-07 -1e-07 0.3000001 0.4763814336453497 1e-07 0 2 1 -4
4 0.0999999 -1.000000000028756e-07 -1e-07 0.3000001 0.2497173224388945 1e-07 0 2 2 -3
1 -1.000000000028756e-07 -1.000000000028756e-07 -1e-07 0.3000001 0.4763814336453497 1e-07 1 3 4 1 4 2 -3
\$EndEntities
\$Nodes
9 190 1 190
0 1 0 1
1
0 0 0
0 2 0 1
2
0.1 0 0
0 3 0 1
3
0.3 0.2 0
0 4 0 1
4
0.3 0.3 0
1 1 0 4
5
6
7
8
0.02000000000000002 0 0
0.04000000000000007 0 0
0.06000000000000008 0 0
0.08000000000000007 0 0
1 2 0 4
9
10
11
12
0.3 0.22 0
0.3 0.24 0
0.3 0.2600000000000001 0
0.3 0.28 0
1 3 0 34
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
0.01605957464567105 0.01180309137910965 0
0.03071752618541084 0.0253154884453655 0
0.04417128399910078 0.04003150743086906 0
0.05661613431529218 0.0556123140212524 0
0.06822049079616135 0.07183008751419329 0
0.07912312645989308 0.08852843782459396 0
0.08943712344795383 0.1055971115351943 0
0.09925499850149465 0.1229561944274134 0
0.1086533091313451 0.1405461651720262 0
0.1176963956684186 0.1583215037111039 0
0.1264393309563879 0.1762464743384211 0
0.134930247471612 0.1942922555620925 0
0.1432122198775771 0.2124349284779746 0
0.1513248457969439 0.2306539942727496 0
0.1593056563070743 0.2489312101134197 0
0.1671914958023691 0.267249621041711 0
0.1750199970366391 0.2855926225365787 0
0.1828313487835908 0.3039429340859255 0
0.1906706587282014 0.3222813173542183 0
0.1985913945697665 0.3405846497580646 0
0.2066609288965863 0.3588228038544409 0
0.2149702701939887 0.3769528336471111 0
0.2236530114099778 0.3949067175619324 0
0.2329274263630457 0.41256101814916 0
0.2432103118394306 0.4296434001974161 0
0.2555297833885212 0.4452870573836619 0
0.2729394211219359 0.4530500479207341 0
0.2854198405876804 0.4383524577409714 0
0.2908727716748887 0.4191972056548036 0
0.2940152340416461 0.399507520698153 0
0.2960445290517715 0.3796687175263951 0
0.2974486675008976 0.3597749756633806 0
0.2984840092269437 0.3398583326125584 0
0.2993025830271605 0.3199314918325419 0
1 4 0 16
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
0.1163729005577897 0.01120840129905042 0
0.1327282389180056 0.02244234005006488 0
0.1489056083611549 0.03393027661042919 0
0.164653950040984 0.04599729734736287 0
0.1796404329120512 0.05899476311588658 0
0.1935319204876485 0.07315383501074177 0
0.2061272577865922 0.08847733934590055 0
0.2174245487270872 0.1047831877538251 0
0.2275715059401495 0.1218306283722394 0
0.2367785348761722 0.1394049843326038 0
0.2452678907478125 0.1573380198155021 0
0.2532705848486707 0.1754941273807766 0
0.2610748370589435 0.193736763607996 0
0.2692178547236704 0.2118288246013934 0
0.2799656870584902 0.2283555559253524 0
0.2935704717799487 0.2187577945522837 0
2 1 0 128
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
0.1949647137023976 0.2873924464840887 0
0.2500849743832158 0.2113860258126266 0
0.2246203816927326 0.1531079440651703 0
0.1976551372003312 0.1064136491457644 0
0.2819915949830282 0.3088240390404488 0
0.1630466926628014 0.2145187274654871 0
0.1610920372499164 0.06494578280807206 0
0.2829797860977996 0.270583618866641 0
0.1461997597002869 0.1795754364756109 0
0.1778416372356152 0.2507436620192215 0
0.1317858844928562 0.1403944846362967 0
0.2807066659979525 0.3488769247886381 0
0.2104099571256749 0.3264419331053582 0
0.1310177067096543 0.04303158259055856 0
0.1093793706670115 0.1057741786791495 0
0.07629739550076865 0.05361027738056431 0
0.2264867913521288 0.3591841032160146 0
0.2782965703200137 0.3878852930667243 0
0.1001267065275372 0.0214085126615364 0
0.07000000000000008 0.01732050807568875 0
0.2417400394284611 0.3962395978689655 0
0.2719695059187498 0.4240453729044958 0
0.2606712198443385 0.228316070193212 0
0.2408288993667751 0.2284848017379404 0
0.2298189458970156 0.2121502900650824 0
0.2511530191307514 0.2447285953702283 0
0.2315016105920924 0.2450702656748411 0
0.2414031576344645 0.2616125606851791 0
0.2217297821240604 0.2617603388757486 0
0.231412231631706 0.2786492918175246 0
0.2507625424359485 0.2786802132158412 0
0.2410667915966609 0.2955680783754726 0
0.2127702069322346 0.2461858540436856 0
0.2228019922748896 0.2951101058456012 0
0.2026242010564789 0.263129247064173 0
0.2633185962632362 0.2976869319987773 0
0.2490069531397094 0.3106885767157005 0
0.2380478284833808 0.1921209581685009 0
0.2172609967636624 0.1964125823620185 0
0.2105841819610411 0.213667115075964 0
0.1984166066693181 0.1982199473280235 0
0.2048061043143432 0.1811350391712491 0
0.1863575872712006 0.1824151415276672 0
0.1937522563372162 0.1646276558517861 0
0.1751561097308549 0.1659868193931924 0
0.1827048639847128 0.1480719515983609 0
0.163527070690436 0.1498577615673423 0
0.1710594454557453 0.1319449998215243 0
0.1933148761200039 0.2161989300566808 0
0.2031894590895256 0.1457609239521086 0
0.1512256870920589 0.1342143599871261 0
0.1592025891276359 0.1159761094205793 0
0.1394322575637569 0.1186055120785693 0
0.1477301490054007 0.1008728891296324 0
0.16716010055073 0.09776197857147706 0
0.2705043813747295 0.2430774805124383 0
0.264257483077097 0.2617292437771941 0
0.1662893148286986 0.1819117314811819 0
0.08791070613985308 0.06995403559351349 0
0.09528857325887263 0.05242496354395079 0
0.1077546872675594 0.06904089583042999 0
0.09893198512275833 0.08786340811265089 0
0.1187890107232253 0.08732012976261749 0
0.1288117270575981 0.06545999025064461 0
0.0843569960305714 0.03323816219263534 0
0.1195796049714526 0.1228774858700743 0
0.06381319458660575 0.03698265952835578 0
0.1380731367051423 0.1615892846814784 0
0.2321936825939065 0.1706620430266402 0
0.2830120920086975 0.2894239197552412 0
0.1696880805436038 0.2333091008742744 0
0.2019580197782724 0.3065323082721415 0
0.04813622415306446 0.0184234596611999 0
0.1788280302340748 0.1134269575722573 0
0.1907514275229403 0.1297713396821061 0
0.1767308718735076 0.201374986158574 0
0.2123105775168673 0.2789526693309205 0
0.2819755576826329 0.3290670210436079 0
0.2635589073920008 0.3384858192984648 0
0.2630038069714074 0.3579780918757836 0
0.2445110812215794 0.3485551741126042 0
0.2459830878649457 0.3673314779650369 0
0.2266943039283714 0.339686512685953 0
0.1158639059466707 0.03093169956414502 0
0.2323153002018318 0.3114759628075434 0
0.2795016325134907 0.3683914561933606 0
0.2626836015013364 0.3765499154013912 0
0.2088211498787052 0.1233509517138604 0
0.1861024137346786 0.2332246805442022 0
0.146630729039449 0.05339665379128253 0
0.2218542514494845 0.2292410093086949 0
0.2551843456681163 0.4134154716157374 0
0.241994580952327 0.3290255727051282 0
0.03098266499682931 0.01110840789713499 0
0.1854580876535818 0.2691706616267356 0
0.1745088327319654 0.07779427457423105 0
0.1555430600106742 0.0830042148715435 0
0.18502134915497 0.09287504969021122 0
0.2753445836402794 0.4067060458275345 0
0.1550681844024139 0.1973513442701535 0
0.1131701986118337 0.04958926058276567 0
0.2688660999765558 0.2796207855227389 0
0.220196206295457 0.1366910864871965 0
0.2601394296130177 0.393720681010417 0
0.1278941521122962 0.1042488679005776 0
0.1376532326817209 0.08488559249989719 0
0.2663834251437626 0.3188447728769553 0
0.2326249113694105 0.3794912978008423 0
0.28467803294503 0.2534615745888956 0
0.1567827226954492 0.1650668898190715 0
0.2119274016186405 0.1648864587486304 0
0.2030494931479687 0.2313789074616253 0
0.1459461572078717 0.07033844684428797 0
0.2914330963462375 0.229555301550553 0
0.100668984559669 0.03687344337662693 0
0.2688919852508933 0.4387236924719528 0
0.2165202743604818 0.3110091585669557 0
0.1481953886746851 0.1491887350472818 0
0.2128047085097481 0.3452001935990733 0
0.08577926330076546 0.01430437695889267 0
0.2581823019731424 0.4292645712284602 0
0.2848342633585693 0.2393966038753567 0
0.1943638700735054 0.2491960793496311 0
0.2483351348595672 0.3823412135158502 0
0.2207983548391477 0.1815016210042254 0
0.2252616414480962 0.3244128170202808 0
0.1771401505218475 0.219591825574952 0
0.2077352925594773 0.2933005365539986 0
\$EndNodes
\$Elements
3 326 1 326
1 1 1 5
1 1 5
2 5 6
3 6 7
4 7 8
5 8 2
1 2 1 5
6 3 9
7 9 10
8 10 11
9 11 12
10 12 4
2 1 2 316
11 155 169 141
12 76 152 126
13 13 1 5
14 97 139 63
15 99 169 155
16 63 157 97
17 126 163 76
18 98 132 67
19 73 115 113
20 67 169 98
21 135 156 6
22 119 171 70
23 154 161 84
24 68 162 138
25 10 176 9
26 73 128 115
27 14 156 135
28 59 100 58
29 118 171 119
30 138 162 120
31 154 166 161
32 7 135 6
33 79 145 143
34 123 126 125
35 70 164 119
36 66 137 136
37 112 165 65
38 65 173 112
39 100 131 58
40 127 129 82
41 137 150 112
42 123 125 124
43 143 144 79
44 66 150 137
45 98 169 99
46 85 86 64
47 90 119 93
48 19 124 77
49 143 155 141
50 64 100 59
51 131 173 65
52 82 135 7
53 150 165 112
54 123 163 126
55 152 175 126
56 131 187 173
57 53 160 52
58 6 156 5
59 98 164 132
60 13 156 14
61 124 125 77
62 18 124 19
63 113 180 73
64 119 164 93
65 60 85 64
66 117 159 158
67 118 119 88
68 82 182 127
69 88 119 90
70 129 135 82
71 78 129 127
72 167 168 116
73 95 185 174
74 85 118 88
75 12 70 11
76 85 88 86
77 121 124 18
78 61 85 60
79 174 185 151
80 12 132 70
81 58 131 57
82 57 131 65
83 136 160 66
84 15 135 129
85 20 128 21
86 61 118 85
87 21 128 73
88 117 160 136
89 77 128 20
90 16 129 78
91 67 132 4
92 23 130 71
93 120 162 71
94 17 121 18
95 66 160 53
96 86 87 64
97 158 160 117
98 62 176 61
99 29 63 30
100 27 133 72
101 63 134 30
102 4 132 12
103 31 134 75
104 22 130 23
105 88 89 86
106 40 84 41
107 26 133 27
108 90 93 92
109 108 137 112
110 46 67 4
111 60 64 59
112 43 148 44
113 93 94 92
114 110 137 108
115 78 122 121
116 30 134 31
117 136 137 110
118 122 123 121
119 87 102 101
120 116 117 114
121 114 136 110
122 78 127 122
123 78 121 17
124 8 82 7
125 27 72 28
126 68 133 26
127 142 143 141
128 117 136 114
129 16 78 17
130 73 130 22
131 108 112 106
132 19 77 20
133 21 73 22
134 25 68 26
135 23 71 24
136 31 75 32
137 115 116 114
138 113 114 110
139 35 83 36
140 42 80 43
141 44 74 45
142 33 79 34
143 109 110 108
144 74 142 141
145 107 108 106
146 123 124 121
147 113 115 114
148 15 129 16
149 57 65 56
150 54 66 53
151 109 113 110
152 14 135 15
153 51 69 50
154 107 109 108
155 49 76 48
156 105 106 104
157 125 168 167
158 47 81 2
159 105 107 106
160 87 100 64
161 102 111 103
162 105 120 107
163 102 103 101
164 111 138 103
165 87 101 100
166 103 138 105
167 103 105 104
168 91 139 97
169 88 90 89
170 103 104 101
171 84 183 154
172 130 172 71
173 93 98 94
174 145 181 75
175 105 138 120
176 74 141 140
177 94 96 92
178 96 139 92
179 142 144 143
180 92 139 91
181 46 140 67
182 48 146 47
183 91 95 89
184 90 91 89
185 45 140 46
186 74 140 45
187 76 146 48
188 90 92 91
189 91 97 95
190 98 99 94
191 47 146 81
192 80 148 43
193 71 172 120
194 94 147 96
195 52 160 158
196 44 148 74
197 29 157 63
198 99 147 94
199 72 185 157
200 55 150 54
201 54 150 66
202 74 148 142
203 72 157 28
204 115 167 116
205 145 155 143
206 127 182 81
207 148 149 142
208 71 162 24
209 128 167 115
210 80 149 148
211 133 151 72
212 77 167 128
213 25 162 68
214 75 188 145
215 86 153 87
216 89 153 86
217 142 149 144
218 99 155 147
219 11 171 10
220 158 159 69
221 50 152 49
222 69 152 50
223 157 185 97
224 76 163 146
225 130 180 172
226 49 152 76
227 9 176 62
228 79 181 145
229 87 153 102
230 36 154 37
231 83 154 36
232 95 153 89
233 79 170 34
234 141 169 140
235 41 161 42
236 144 170 79
237 65 165 56
238 132 164 70
239 102 174 111
240 126 168 125
241 35 170 83
242 112 173 106
243 28 157 29
244 120 172 107
245 84 161 41
246 116 159 117
247 52 158 51
248 51 158 69
249 100 187 131
250 2 182 8
251 81 182 2
252 73 180 130
253 42 161 80
254 24 162 25
255 101 187 100
256 125 167 77
257 111 174 151
258 122 163 123
259 75 181 32
260 163 177 146
261 33 181 79
262 151 189 111
263 133 189 151
264 56 165 55
265 83 166 154
266 80 166 149
267 70 171 11
268 93 164 98
269 138 189 68
270 55 165 150
271 171 184 10
272 34 170 35
273 140 169 67
274 81 177 127
275 166 186 149
276 122 177 163
277 97 185 95
278 161 166 80
279 111 189 138
280 116 168 159
281 83 186 166
282 107 172 109
283 106 173 104
284 61 184 118
285 37 183 38
286 5 156 13
287 39 178 40
288 40 178 84
289 69 175 152
290 155 188 147
291 146 177 81
292 172 180 109
293 173 187 104
294 179 190 96
295 127 177 122
296 134 179 75
297 153 174 102
298 95 174 153
299 147 179 96
300 63 190 134
301 38 178 39
302 134 190 179
303 8 182 82
304 159 175 69
305 32 181 33
306 109 180 113
307 96 190 139
308 139 190 63
309 104 187 101
310 154 183 37
311 68 189 133
312 145 188 155
313 170 186 83
314 168 175 159
315 126 175 168
316 149 186 144
317 151 185 72
318 144 186 170
319 118 184 171
320 10 184 176
321 38 183 178
322 147 188 179
323 179 188 75
324 178 183 84
325 176 184 61
326 62 3 9
\$EndElements
``````

Here’s the FEniCSx code:

``````def thermal_problem(meshfile):
current_curves=(1,2)
voltages_curves=(1,2)
inj_current_curve, out_current_curve = current_curves

# Define FE function space
deg = 3
el = FiniteElement("CG", mesh.ufl_cell(), deg)
V = FunctionSpace(mesh, el)

u = TestFunction(V)
temp = Function(V)

κ = 1.8
ΔT = 30
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.

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}ksp_type"] = "preonly"
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}''')

# Compute fluxes on boundaries
n = FacetNormal(mesh)
down_heat_flux = form(dot(-κ * grad(temp), n)*ds(out_current_curve))
Q_out = assemble_scalar(down_heat_flux)

top_heat_flux = form(dot(-κ * grad(temp), n)*ds(inj_current_curve))
Q_in = assemble_scalar(top_heat_flux)

print(f'''------- Post processing --------
Q_in: {Q_in}
Q_out: {Q_out}''')
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)
with VTXWriter(MPI.COMM_WORLD, "results/heat_flux.bp", [q], engine="BP4") as vtx:
vtx.write(0.0)
return Q_in, Q_out

t_pb = thermal_problem("meshes/mesh.msh")
``````

It gives me the output

``````Info    : Reading 'meshes/mesh.msh'...
Info    : 17 entities
Info    : 190 nodes
Info    : 326 elements
------- Pre-processing --------
Length of the side where heat enters: 0.1
Length of the side where heat leaves the material: 0.09999999999999998

------- Processing --------
Number of interations: 1
------- Post processing --------
Q_in: 13.764246963696966
Q_out: -13.887692194019193
``````

Which is not acceptable (Q_in should be equal to Q_out, at least much closer), I picked a FE space of degree 3, and even with degree 6 I still don’t get matching fluxes.

In Paraview, I see:

And Gmsh shows the mesh:

You are moving a Discontinuous function into a continuous space, which will cause approximation issues, see
https://jsdokken.com/FEniCS23-tutorial/src/approximations.html

It is still not clear what you want to use `q` for.
Note that `q` doesn’t depend on the FacetNormal.

That it is in a test file means that it is tested for correctness. It is usable with the main branch of dolfinx, ffcx, Basix and ufl.

1 Like

I don’t understand why that would be the case? My code is:

``````    # Define FE function space
deg = 3
el = FiniteElement("CG", mesh.ufl_cell(), deg)
V = FunctionSpace(mesh, el)
``````

And q involves grad(temp), I assume that internally FEniCSs converts it as a vector in a CG space of degree 2. But it’s DG somewhat, and which degree?

I am mostly interested in getting Q_in and Q_out accurately, as their difference is a major “result” of the simulations in my case. And this process involves facetnormal.

I build q because I want to visualize how the heat flux (as close to FEniCSx internal as possible) looks like, to try to understand why my Q_in - Q_out is faulty, wrong, and still not converged when I use a FE degree of 6 (on a denser mesh than my example here).

Edit: About your link of Sorbonne, the take away is that using a DG space where the discontinuity is aligned with the mesh “fixes” the problem?

Temp only has continuous higher order derivatives within a single element. At a vertex between multiple cells, the gradient will be different when viewed at a single point if viewed from either cell.

`Q_in` and `Q_out` uses integration over facets, which are well defined for standard quadrature schemes.

How does Q_in and Q_out behave on a finer grid?

For what I consider a high density mesh, I get

``````Info    : Reading 'meshes/mesh.msh'...
Info    : 17 entities
Info    : 27963 nodes
Info    : 55335 elements
------- Pre-processing --------
Length of the side where heat enters: 0.1
Length of the side where heat leaves the material: 0.09999999999999998

------- Processing --------
Number of interations: 1
------- Post processing --------
Q_in: 16.830789755954694
Q_out: -16.917597520767256
(16.830789755954694, -16.917597520767256)
``````

OK about your comment about the issue to pass from a DG to CG. I now built q using a DG space of degree 2, for the visualization in Paraview. If I scale by magnitude, and I zoom into a corner, this is what I see:

Hmm these big magnitude arrows have the same origin…

That is not suprising, as a DG space has duplicate dofs at each vertex. This usually means that the mesh should be further refined in that area, as the variation between two cells are really large.

Just a side note, the coarse and fine mesh you have presented look very different.

Thanks for all dokken. Yes, everytime I generate a mesh, it changes (I use randomness to create 2 Bezier curves between the 2 fixed curves where heat enters. I will modify it with an algorithm as to maximize Q_in-Q_out in a more complicated problem than the pure thermal one. The algorithm will be guided by Q_in-Q_out.
Ok about the mesh refinement, thanks for your input, I think I can easily improve the refinement locally there using Gmsh.

I am still not entirely satisfied with the result of Q_in-Qout, the result differ by 1 in 168. Any other computation I perform (like the electrical resistance) is much more accurate than 1 part per 168.

But I will try with a smarter mesh. If you have any other suggestion, please let me know.

Hello @dokken I think I have learned a few things and I felt like sharing with you (and all others who are interested). I noticed, after many tries to play with a mesh with a high resolution near the 2 sides where the heat enters/leaves (i.e. where I compute the fluxes), that Q_in changes relatively very little when either the mesh is refined or the CG space of “temp” is increased, whereas Q_out changes a lot.

This makes me believe that the fluxes are reasonably well computed, on both sides. The fact that Q_in-Q_out hasn’t converged is an indication that the solution itself hasn’t converged properly. So, a fix would be to refine the whole mesh, I think, and possibly increase the degree of the FE space of `temp`.

dokken, do you have any idea if it is better to use odd, as opposed to even, order for a FE space? If I choose an odd order, then the flux will be defined in an even DG space. And since I take derivatives of this quantity with respect to position, I guess I end up with another DG space of 1 degree less, for that quantity.

Hello @dokken, I have a new question.I have just read that the fluxes are supposed to be exact only at quadrature points and that if I want for example to evaluate the total flux accurately through a surface, I need to employ a method dealing with virtual work of internal forces. Source: Computing consistent reaction forces — Computational Mechanics Numerical Tours with FEniCSx

I wanted to know your opinion, i.e. is it really true? Could this ‘‘solve’’ my problem regarding Qin and Qout?
If so, I’ll open a new thread as my code is more complex (mixed elements) than the one here, and I have some doubts on how to implement.this.

I trust Jeremy. If he states that it is the only way to accurately evaluate the fluxes that is probably true. It could probably solve your issue.

1 Like

Thanks for your opinion dokken.I tried a few things without success so far, I will post a MWE when I have more time.

I do wonder though whether I can apply this strategy to a thermal flux evaluation (there is no '‘virtual work’'of thermal flux in a solid, though I am unsure if there exist an equivalent concept). The virtual force in this case would be grad(T).