Computing fluxes accurately

mathematics
1

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?

2

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.

3

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:
    # Read in mesh
    msh = meshio.read('the_mesh.msh')
    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 = xdmf.read_mesh(name="Grid")
    ct = xdmf.read_meshtags(mesh, name="Grid")

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

# 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 %).
Figure_1
vs
Figure_1

4

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

P1
P13136×1044 290 KB

P2
P23136×1044 235 KB

P3
P33136×1044 225 KB

5

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.

6

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).

7

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?

8

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?

9

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
10

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
    reading_voltage_curve_0, reading_voltage_curve_1 = voltages_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.
    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}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
Info    : Done reading 'meshes/mesh.msh'
Info    : Reading 'meshes/mesh.msh.opt'...
 ------- 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:

heat_flux
heat_flux1367×789 73.7 KB

And Gmsh shows the mesh:
whale
whale800×915 89.6 KB

11

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
12

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?

13

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?

14

For what I consider a high density mesh, I get

Info    : Reading 'meshes/mesh.msh'...
Info    : 17 entities
Info    : 27963 nodes
Info    : 55335 elements
Info    : Done reading 'meshes/mesh.msh'
Info    : Reading 'meshes/mesh.msh.opt'...
 ------- 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)

whale2
whale2800×915 158 KB

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:
heat_flux_2
heat_flux_21367×789 100 KB

Hmm these big magnitude arrows have the same origin…

15

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.

16

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.

17

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.

18

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.

19

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
20

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).

I wonder if @bleyerj has an idea about this.