Coverage for pyxfem/core/finite_volume_2d.py: 0%
241 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-10-24 07:17 +0000
1#!/usr/bin/env python3
2# -*- coding: utf-8 -*-
3# pylint: disable=C0103, W1401, C0302
4r"""
5Created on Thu Mar 7 11:17:17 2024
6@author: nchevaug
72D Explicit Lagrangian Finite Volume Solver (GLACE Scheme).
9Philosophy & Architecture:
10This module implements the Godunov-type Lagrangian Scheme for Compressible
11Euler equations (GLACE). It computes the evolution of a fluid on an unstructured
122D triangular mesh moving with the fluid velocity (Lagrangian framework).
14The solver features a dual-execution architecture:
151. Pure NumPy (Fallback): Fully vectorized array operations.
162. Numba JIT (HPC Mode): Aggressively optimized loop-level parallelism
17 (`prange`), avoiding heap allocations and leveraging hardware registers
18 for the nodal Riemann solver.
19"""
20from abc import abstractmethod
21from typing import Tuple, NamedTuple, Sequence, List, Union, assert_never, Optional, overload
22#from typing import Any, Callable, Dict
23#from collections import defaultdict
24from llvmlite.ir import FastMathFlags
25import numpy as np
26import numba
27#from numba import guvectorize, float32, float64, prange
29from pyxfem.core.config import nbfloat
30from pyxfem.core.nbkalgebra2d import Vec2, Sy22
31import numba
32import numba.core.types as nbct
33import numba.typed as nbty
34# pylint: disable=ungrouped-imports
35from pyxfem.core.nbkalgebra2d import \
36 vec2_t, sy22_t, vec2_zero, sy22_zero,\
37 set_vec2,\
38 vec2_scal, vec2_square, vec2_norm, vec2_dot, vec2_add, vec2_add3, vec2_sub, \
39 sy22_add, sy22_mv, sy22_smv, sy22_r, sy22_r2, sy22_trace
42from pyxfem.core.graph import T2V, extract_txy, upward_adjacencies, which_in_upward
43from pyxfem.core.config import USE_NUMBA, FASTMATH, IndexType, npfloat, fzero, fone, fhalf,\
44 FloatArray, IndexArray
45from pyxfem.core.smallsize_tensor import max_eigval22sym, gemm22, inv, det, sysmv22, sysmv22_guvec, vec_norm2
46import pyxfem.core.graph as graph
48def t_tangentmap(tris_coord, out: Optional[FloatArray] = None):
49 ''' return the "tangent map"
50 The F tensor between reference element and each of the elment in tris_coord
51 '''
52 shape = tris_coord.shape[:-2] +(2,2)
53 out = np.empty(shape, npfloat) if out is None else out
54 v0 = tris_coord[..., 0, :]
55 out[... , 0] = tris_coord[..., 1, :] - v0
56 out[..., 1] = tris_coord[..., 2, :] - v0
57 #outold = tris_coord[...,1:,:] -tris_coord[..., 0, np.newaxis, :]
58 #print(np.linalg.norm(out - outold.swapaxes(-1,-2)))
59 return out
61 return np.subtract(tris_coord[...,1:,:], tris_coord[..., 0, np.newaxis, :], out = out)
64class GELEES_MeshStruct(NamedTuple):
65 """
66 A structured container for the 2D Lagrangian mesh topology, geometry,
67 and inverse connectivities.
69 This NamedTuple holds all the static and dynamically updated arrays
70 required for the finite volume solver, designed to be easily passed
71 into Numba JIT-compiled functions.
73 Attributes
74 ----------
75 ncells : int
76 The total number of triangular cells in the mesh.
77 nverts : int
78 The total number of vertices (nodes) in the mesh.
79 t2v : IndexArray
80 2D array of shape `(ncells, 3)`. The cell-to-vertex connectivity matrix
81 mapping each triangle to its 3 defining vertex indices.
82 v2t_index : IndexArray
83 1D array of shape `(nverts + 1,)`. The index pointer array for the CSR-like
84 inverse connectivity (vertices-to-cells).
85 v2t_tris : IndexArray
86 1D array containing the flattened cell indices surrounding each vertex.
87 Sliced using `v2t_index`.
88 v2t_v : IndexArray
89 1D array containing the local vertex index (0, 1, or 2) within the
90 corresponding cell in `v2t_tris`. Sliced using `v2t_index`.
91 xy : FloatArray
92 2D array of shape `(nverts, 2)`. The physical (x, y) coordinates of
93 each vertex. This array is updated at each time step.
94 S : FloatArray
95 1D array of shape `(ncells,)`. The surface area of each triangular cell.
96 C : FloatArray
97 3D array of shape `(ncells, 3, 2)`. The corner normal vectors (C_jr)
98 for each vertex of each cell.
99 nC : FloatArray
100 2D array of shape `(ncells, 3)`. The Euclidean norms (magnitudes) of
101 the corner normal vectors.
102 N : FloatArray
103 3D array of shape `(ncells, 3, 2)`. The normalized unit direction
104 vectors of the corner normals.
105 """
106 ncells: int
107 nverts: int
108 t2v: T2V
109 v2t_index: IndexArray
110 v2t_tris: IndexArray
111 v2t_v: IndexArray
112 xy: FloatArray
113 cell_meas: FloatArray
114 C: FloatArray
115 nC: FloatArray
116 N: FloatArray
118def allocate_Cjr(nc :int):
119 r"""
120 Allocates memory for the geometric properties of the mesh cells.
122 This function initializes empty NumPy arrays required to store the
123 surface areas and the corner normal vectors (and their properties)
124 for a 2D triangular mesh.
126 Parameters
127 ----------
128 nc : int
129 The total number of triangular cells in the mesh.
131 Returns
132 -------
133 S : ndarray
134 1D array of shape `(nc,)` and type `float64`.
135 Allocated to store the surface area of each triangular cell.
136 C : ndarray
137 3D array of shape `(nc, 3, 2)` and type `float64`.
138 Allocated to store the corner normal vectors (C_jr) for the 3 vertices
139 of each cell in 2D space (x, y).
140 normC : ndarray
141 2D array of shape `(nc, 3)` and type `float64`.
142 Allocated to store the Euclidean norms (magnitudes) of the C_jr vectors.
143 """
144 cell_measure = np.empty( (nc), dtype = npfloat)
145 C = np.empty( (nc, 3, 2), dtype = npfloat)
146 normC = np.empty( (nc, 3), dtype = npfloat)
147 N = np.empty( (nc, 3, 2), dtype = npfloat)
148 return cell_measure, C, normC, N
150def setup_MeshStruct(xy, tris):
151 v2t = upward_adjacencies(tris)
152 v2t_v = which_in_upward( v2t, tris)
153 ncells = tris.shape[0]
154 nverts = xy.shape[0]
155 cell_meas, C, normC, N = allocate_Cjr(ncells)
156 mesh = GELEES_MeshStruct(
157 ncells =ncells, nverts=nverts,
158 t2v=tris, v2t_index=v2t.index,
159 v2t_tris= v2t.data, v2t_v=v2t_v,
160 xy=xy,
161 cell_meas=cell_meas,
162 C=C, nC=normC, N=N)
163 if USE_NUMBA:
164 computeCjr_triangles_numba(mesh)
165 else:
166 computeCjr_triangles(mesh)
167 return mesh
169def update_MeshStruct(xy, mesh : GELEES_MeshStruct):
170 mesh.xy[:] = xy
171 computeCjr_triangles_numba(mesh)
173def computeCjr_triangles(mesh : GELEES_MeshStruct):
174 xy, tris = mesh.xy, mesh.t2v
175 txy = extract_txy(xy, tris)
176 F1 = t_tangentmap(txy)
177 computeCjr_triangles_(F1, mesh.cell_meas, mesh.C, mesh.nC, mesh.N)
179def computeCjr_triangles_(F, outMeas, outC, outnormC, outN):
180 det(F, outMeas)
181 outMeas*=0.5
182 F11 = F[...,0,0]
183 F12 = F[...,0,1]
184 F21 = F[...,1,0]
185 F22 = F[...,1,1]
186 #outC[:] = fhalf*np.stack([-F22+F12, F21-F11, F22, -F21, -F12, F11], axis=-1).reshape((-1, 3, 2))
187 np.subtract(F12, F22, outC[..., 0,0])
188 np.subtract(F21, F11, outC[..., 0,1])
189 np.copyto(outC[...,1,0], F22)
190 np.multiply(-fone, F21, outC[...,1,1])
191 np.multiply(-fone, F12, outC[...,2,0])
192 np.copyto(outC[...,2,1], F11)
193 np.multiply(fhalf, outC, outC)
194 vec_norm2(outC, outnormC)
195 np.divide(outC, outnormC[..., None], out = outN)
197@numba.njit(inline='always', fastmath = FASTMATH)
198def get_vec2(index, array):
199 return (array[index, 0], array[index, 1])
201@numba.njit(cache=True, fastmath = FASTMATH)
202def get_trixy(elemid, mesh):
203 tri = mesh.t2v[elemid]
204 n0, n1, n2 = tri[0], tri[1], tri[2]
205 xy = mesh.xy
206 xy0 = get_vec2(n0, xy)
207 xy1 = get_vec2(n1, xy)
208 xy2 = get_vec2(n2, xy)
209 return xy0, xy1, xy2
211@numba.njit(cache = True, fastmath = FASTMATH)
212def Cjr_triangle(xy0: Vec2, xy1: Vec2, xy2: Vec2):
213 """
214 Computes the surface area and corner normal vectors (C_jr) for a single
215 2D triangle.
217 This function is optimized for Numba JIT compilation. It performs
218 computations using primitive types and returns tuples to ensure
219 zero-allocation overhead when called inside high-performance loops
220 (like `prange`).
222 Parameters
223 ----------
224 xy0 : tuple of float or array-like
225 The (x, y) coordinates of the 0-th vertex of the triangle.
226 xy1 : tuple of float or array-like
227 The (x, y) coordinates of the 1-st vertex of the triangle.
228 xy2 : tuple of float or array-like
229 The (x, y) coordinates of the 2computeCjr_triangles_numba(mesh)
231def computeCjr_triangles(F, outMeas, outC, outnormC, outN):
232 -nd vertex of the triangle.
234 Returns
235 -------
236 S : float
237 The surface area of the triangle.
238 Cj0 : tuple of float
239 A tuple `(x, y)` representing the corner normal vector
240 associated with vertex 0.
241 Cj1 : tuple of float
242 A tuple `(x, y)` representing the corner normal vector
243 associated with vertex 1.
244 Cj2 : tuple of float
245 A tuple `(x, y)` representing the corner normal vector
246 associated with vertex 2.
247 """
248 F0 = (xy1[0] - xy0[0] , xy1[1] - xy0[1])
249 F1 = (xy2[0] - xy0[0] , xy2[1] - xy0[1])
250 S = npfloat(0.5)*F0[0]*F1[1] - npfloat(0.5)*F1[0]*F0[1]
251 Cj0 = (npfloat(0.5)*(F0[1]-F1[1]), npfloat(0.5)*(F1[0]-F0[0]))
252 Cj1 = (npfloat(0.5)*F1[1], -npfloat(0.5)*F1[0])
253 Cj2 = (-npfloat(0.5)*F0[1], npfloat(0.5)*F0[0])
254 return S, Cj0, Cj1, Cj2
257@numba.njit(parallel=True, cache = True, fastmath = FASTMATH)
258def computeCjr_triangles_numba(mesh):
259 """
260 Computes and updates the geometric properties of all triangular cells
261 in the mesh using Numba parallelization.
263 This function iterates over all cells concurrently, calculates the
264 surface area and the corner normal vectors (C_jr) for each triangle
265 based on its current vertex coordincomputeCjr_triangles_numba(mesh)
267 Parameters
268 ----------
269 mesh : NamedTuple (e.g., GELEES_MeshStruct)
270 A structured object containing the mesh data. It must provide:
271 - `mesh.ncells` (int): The total number of cells.
272 - `mesh.S` (ndarray): 1D array of shape `(ncells,)` to store updated areas.
273 - `mesh.C` (ndarray): 3D array of shape `(ncells, 3, 2)` to store
274 updated corner normal vectors.
275 - Connectivity and coordinate data required by `get_trixy`.
277 Returns
278 -------
279 None
280 The function does not return any value. The `mesh.S` and `mesh.C`
281 arrays are modified in-place for maximum performance.
282 """
283 N = mesh.ncells
284 outS = mesh.cell_meas
285 outC = mesh.C
286 outnC = mesh.nC
287 outN = mesh.N
288 for j in numba.prange(N): # type ignore # pylint: disable=E1133
289 xy0, xy1, xy2 = get_trixy(j, mesh)
290 Sj, Cj0, Cj1, Cj2 = Cjr_triangle( xy0, xy1, xy2)
291 nCj0, nCj1, nCj2 = \
292 vec2_norm(Cj0), vec2_norm(Cj1), vec2_norm(Cj2)
293 outS[j] = Sj
294 outnC[j, 0] = nCj0
295 outnC[j, 1] = nCj1
296 outnC[j, 2] = nCj2
297 for k in range(2):
298 outC[j, 0, k] = Cj0[k]
299 outC[j, 1, k] = Cj1[k]
300 outC[j, 2, k] = Cj2[k]
301 outN[j, 0, k] = Cj0[k]/nCj0
302 outN[j, 1, k] = Cj1[k]/nCj1
303 outN[j, 2, k] = Cj2[k]/nCj2
305class ConservativeLawManager():
306 def __init__(self, nbdata, U_shapes, U_indexes, W_shapes, W_indexes):
307 self.U_indexes = U_indexes
308 self.U_shapes = U_shapes
309 sizeU = U_shapes["size"]
310 self.W_indexes = W_indexes
311 self.W_shapes = W_shapes
312 sizeW = W_shapes["size"]
313 self.U = np.zeros((nbdata, sizeU), dtype=npfloat)
314 self.W = np.zeros((nbdata, sizeW), dtype=npfloat)
316 @abstractmethod
317 def conservative_to_primal(self, cell_measure):
318 pass
320 @abstractmethod
321 def acoustic_impedance_tensor(self, n, outZ = None):
322 pass
324 @overload
325 def get_primal_variable(self, names: str) -> FloatArray: ...
326 @overload
327 def get_primal_variable(self, names: Sequence[str]) -> List[FloatArray]: ...
328 def get_primal_variable(self, names: Union[str, Sequence[str]]) -> FloatArray | List[FloatArray]:
329 """ return primal variable corresponding to it's name or a list of names; avalable variables are listed in W_indices """
330 def get(name:str):
331 s, e = self.W_indexes[name]
332 return self.W[:, s:e].reshape((-1,)+self.W_shapes[name])
333 match names :
334 case str() as name:
335 return get(name)
336 case [*sequ_names] :
337 res = []
338 for name in sequ_names:
339 res+=[get(name)]
340 return res
341 case _ as unreachable_var:
342 assert_never(unreachable_var)
343 @overload
344 def get_conservative_variable(self, names: str) -> FloatArray: ...
345 @overload
346 def get_conservative_variable(self, names: Sequence[str]) -> List[FloatArray]: ...
348 def get_conservative_variable(self, names: Union[str, Sequence[str]]) -> FloatArray | List[FloatArray]:
349 """ return primal variable corresponding to it's name or a list of names; avalable variables are listed in W_indices """
350 def get(name:str):
351 s, e = self.U_indexes[name]
352 return self.U[:, s:e].reshape((-1,)+self.U_shapes[name])
353 match names :
354 case str() as name:
355 return get(name)
356 case [*sequ_names] :
357 res = []
358 for name in sequ_names:
359 res+=[get(name)]
360 return res
361 case _ as unreachable_var:
362 assert_never(unreachable_var)
364class GlaceExplicitLagrangianBaseSolver:
365 def __init__(self, conservation_law : ConservativeLawManager, mesh: GELEES_MeshStruct):
366 self.conservation_law = conservation_law
367 self.mesh = mesh
368 nverts = mesh.nverts
369 ncells = mesh.ncells
370 self.vA = np.zeros((nverts, 2, 2), dtype=npfloat)
371 self.vb = np.zeros((nverts, 2), dtype=npfloat)
372 self.vu = np.zeros((nverts, 2), dtype=npfloat)
374 self.cvp = np.zeros((ncells, 3), dtype=npfloat)
375 self.cvqflux = np.zeros((ncells, 3, 2), dtype=npfloat)
377 self.cvA = np.zeros((ncells, 3, 2, 2), dtype=npfloat)
378 self.cvb = np.zeros((ncells, 3, 2), dtype=npfloat)
379 self.cvu = np.zeros((ncells, 3, 2), dtype=npfloat)
380 self.qflux = np.zeros((ncells, 2), dtype=npfloat)
381 self.eflux = np.zeros((ncells), dtype=npfloat)
383 def _riemann_solver_assemble(self)-> None:
384 self.vA.fill(fzero)
385 self.vb.fill(fzero)
386 # order = self.mesh.t2v.ravel()
387 # np.add.at(self.vA, order, self.cvA.reshape((-1,2,2)))
388 # np.add.at(self.vb, order, self.cvb.reshape((-1,2)))
389 GlaceExplicitLagrangianBaseSolver._riemann_solver_assemble_numba(self.cvA, self.cvb, self.mesh, self.vA, self.vb)
390 self._set_boundary_conditions()
392 @staticmethod
393 @numba.njit(parallel=True, cache=True, fastmath = FASTMATH)
394 def _riemann_solver_assemble_numba(cvA, cvb, mesh, out_vA, out_vb)-> None:
395 for vid in numba.prange(mesh.nverts): # pylint: disable=E1133
396 s = mesh.v2t_index[vid]
397 e = mesh.v2t_index[vid+1]
398 Ar : Sy22= (out_vA[vid,0,0] ,out_vA[vid,0,1], out_vA[vid,1,1])
399 br : Vec2 = (out_vb[vid,0], out_vb[vid,1])
400 # Cell k connectd to node vid
401 for k in range(s, e):
402 tj = mesh.v2t_tris[k]
403 vr = mesh.v2t_v[k]
404 Ar = sy22_add(Ar, (cvA[tj, vr, 0,0], cvA[tj, vr, 0,1], cvA[tj, vr, 1,1]))
405 br = vec2_add(br, (cvb[tj, vr, 0], cvb[tj, vr, 1]))
406 out_vA[vid, 0, 0] = Ar[0]
407 out_vA[vid, 0, 1] = Ar[1]
408 #out_vA[vid, 1, 0] = Ar[1]
409 out_vA[vid, 1, 1] = Ar[2]
410 set_vec2(vid, out_vb, br)
412 def _riemann_solver(self)-> Tuple[FloatArray, FloatArray]:
413 #self._riemann_solver_setup_0()
414 #self._riemann_solver_setup_1()
415 self._riemann_solver_setup()
416 self._riemann_solver_assemble()
417 sysmv22(self.vA, self.vb, self.vu, parallel =True)
418 graph.extract_txy(self.vu, self.mesh.t2v, self.cvu)
419 self._riemann_solver_cvqflux()
420 return self.vu, self.cvqflux
422 def _riemann_solver_cvqflux(self):
423 if USE_NUMBA:
424 GlaceExplicitLagrangianBaseSolver._riemann_solver_cvqflux_njit(self.cvA, self.cvu, self.cvb, self.cvqflux)
425 else:
426 np.matmul(self.cvA, self.cvu[..., np.newaxis], out=self.cvqflux[..., np.newaxis])
427 np.subtract(self.cvb, self.cvqflux, out=self.cvqflux)
429 @staticmethod
430 @numba.njit(parallel=True, cache = True)
431 def _riemann_solver_cvqflux_njit(cvA, cvu, cvb, out_cvqflux):
432 for j in numba.prange(cvA.shape[0]): # pylint: disable=E1133
433 for r in numba.prange(cvA.shape[1]): # pylint: disable=E1133
434 Ajr : Sy22 = (cvA[j,r, 0, 0 ] , cvA[j,r, 0, 1 ], cvA[j,r, 1, 1 ])
435 ujr : Vec2 = (cvu[j,r, 0], cvu[j,r,1])
436 bjr : Vec2 = (cvb[j,r, 0], cvb[j,r,1])
437 tmpjr = sy22_mv(Ajr, ujr)
438 qjr = vec2_sub(bjr, tmpjr)
439 out_cvqflux[j,r, : ] = qjr
441 @staticmethod
442 @numba.njit(parallel=True, cache=True)
443 def critical_timestep_pnjit(cell_measure, nC, c):
444 min_dt = np.inf # numba does autoamatique reduction on this variable
445 for j in numba.prange(len(cell_measure)): # pylint: disable=E1133
446 cell_measurej = cell_measure[j]
447 nCj0 = nC[j,0]
448 nCj1 = nC[j,1]
449 nCj2 = nC[j,2]
450 cj = c[j]
451 sum_normCj = nCj0 + nCj1 + nCj2
452 dt = cell_measurej / sum_normCj/ cj
453 # # Numba comprend automatiquement qu'il s'agit d'une réduction parallèle
454 min_dt = min(min_dt, dt)
455 return min_dt
457 @overload
458 def get_primal_variable(self, names: str) -> FloatArray: ...
460 @overload
461 def get_primal_variable(self, names: Sequence[str]) -> List[FloatArray]: ...
463 def get_primal_variable(self, names: Union[str, Sequence[str]]) -> FloatArray | List[FloatArray]:
464 """ return primal variable corresponding to it's name or a list of names; avalable variables are listed in W_indices """
465 return self.conservation_law.get_primal_variable(names)
467 @overload
468 def get_conservative_variable(self, names: str) -> FloatArray: ...
470 @overload
471 def get_conservative_variable(self, names: Sequence[str]) -> List[FloatArray]: ...
473 def get_conservative_variable(self, names: Union[str, Sequence[str]]) -> FloatArray | List[FloatArray]:
474 return self.conservation_law.get_conservative_variable(names)