Coverage for pyxfem/fluidmechanics/euler2d_finite_volume.py: 0%
521 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 collections import defaultdict
21from typing import Any, Callable, Dict, Tuple, NamedTuple, Sequence, List, Union, assert_never, Optional, overload
23from numba import guvectorize, float32, float64
25from numba.np.ufunc import parallel
26import numpy as np
28from pyxfem.core.config import (
29 FASTMATH,
30 FloatArray,
31 IndexArray,
32 IndexType,
33 USE_NUMBA,
34 fhalf,
35 fone,
36 npfloat,
37)
38from pyxfem.fluidmechanics.exact_riemann_solver import PerfectGaz
39from pyxfem.core.finite_volume_2d import (
40 ConservativeLawManager,
41 GELEES_MeshStruct,
42 setup_MeshStruct,
43 update_MeshStruct,
44 GlaceExplicitLagrangianBaseSolver,
45 Cjr_triangle
46)
48if USE_NUMBA:
49 from pyxfem.core.config import nbfloat
50 from pyxfem.core.nbkalgebra2d import Vec2, Sy22
51 import numba
52 import numba.core.types as nbct
53 import numba.typed as nbty
54 # pylint: disable=ungrouped-imports
55 from pyxfem.core.nbkalgebra2d import \
56 vec2_t, sy22_t, vec2_zero, sy22_zero,\
57 vec2_scal, vec2_square, vec2_norm, vec2_dot, vec2_add, vec2_add3, vec2_sub, \
58 sy22_add, sy22_mv, sy22_smv, sy22_r, sy22_r2, sy22_trace
60# --- INDICES DES VARIABLES CONSERVATIVES (Tableau U) ---
61DIM = 2
62I_M = 0 # Masse (m)
63I_MU_START = 1
64I_MU_END = I_MU_START + DIM
65I_ME = I_MU_END # Énergie Totale (mE)
66SIZE_U = I_ME+1
68# --- INDICES DES VARIABLES PRIMITIVES (Tableau W, optionnel) ---
69I_RHO = 0 # Densité
70I_U_START = 1 # Vitesse X
71I_U_END = I_U_START + DIM # Vitesse Y
72I_P = I_U_END # Pression
73I_C = I_P+1 # Sound speed
74SIZE_W = I_C+1
76U_indexes = {"m": (I_M, I_M+1), "mu":(I_MU_START, I_MU_END), "mE":(I_ME, I_ME+1) }
77U_shapes = {"size": SIZE_U, "m": (), "mu":(DIM,), "mE":() }
79W_indexes = {"rho": (I_RHO, I_RHO+1), "u": (I_U_START, I_U_END), "p":(I_P, I_P+1), "c":(I_C, I_C+1)}
80W_shapes = {"size":SIZE_W, "rho": (), "u":(DIM,), "p":(), "c":()}
83@guvectorize([
84 (float64, float64, float64[:], float64[:,:]),
85 (float32, float32, float32[:], float32[:,:])
86], '(),(),(m)->(m,m)', nopython=True, cache=True)
87def _acoustic_impedance_tensor_2d_guvec(rho, c, n, out):
88 out[0,0] = rho*c*n[0]*n[0]
89 out[0,1] = rho*c*n[0]*n[1]
90 out[1,0] = rho*c*n[1]*n[0]
91 out[1,1] = rho*c*n[1]*n[1]
93@numba.njit(parallel=True, cache = True)
94def _acoustic_impedance_tensor_2d_njit(rho, c, n, out):
95 N = len(rho)
96 nv = n.shape[-2]
97 for j in numba.prange(N): # type ignore # pylint: disable=E1133
98 rhocj = rho[j] * c[j]
99 for r in range(nv):
100 njr0 = n[j,r,0]
101 njr1 = n[j,r,1]
102 out[j,r,0,0] = rhocj*njr0*njr0
103 out[j,r,0,1] = rhocj*njr0*njr1
104 out[j,r,1,0] = rhocj*njr1*njr0
105 out[j,r,1,1] = rhocj*njr1*njr1
108class ConservationLawManager_Euler(ConservativeLawManager):
109 def __init__(self, gazlaw, nbdata):
110 super().__init__(nbdata, U_shapes, U_indexes, W_shapes, W_indexes)
111 self.gazlaw = gazlaw
112 if USE_NUMBA:
113 self.numba_pressure_c_func, self.numba_soundspeed_c_func, self.fluidlaw_data = self.gazlaw.get_numba_c_funcs()
115 def conservative_to_primal(self, cellmeasure):
116 m, mu, mE = self.get_conservative_variable(["m","mu","mE"])
117 rho, u, p, c = self.get_primal_variable(["rho","u", "p", "c"])
118 if USE_NUMBA :
119 pfunc = self.numba_pressure_c_func
120 cfunc = self.numba_soundspeed_c_func
121 fluidlaw_data =self.fluidlaw_data
122 ConservationLawManager_Euler.conservative_to_primal_numba(m, mu, mE, cellmeasure, pfunc, cfunc, fluidlaw_data, rho, u, p, c)
123 else:
124 np.divide(m, cellmeasure, out = rho)
125 np.divide(mu, m[..., np.newaxis], out = u)
126 u2 = np.square(u)
127 e = mE/m-fhalf*(u2[...,0] + u2[...,1])
128 self.gazlaw.pressure(e, rho, out = p)
129 self.gazlaw.soundspeed(p, rho, out = c)
131 if USE_NUMBA :
132 @staticmethod
133 @numba.njit(parallel=True, cache = True)
134 def conservative_to_primal_numba(m, mu, mE, cell_meas, pfunc, cfunc, fluidlaw_data, rho, u, p, c):
135 N = len(m)
136 for i in numba.prange(N): # type ignore # pylint: disable=E1133
137 rho[i] = m[i]/cell_meas[i]
138 u[i, 0] = mu[i,0]/m[i]
139 u[i, 1] = mu[i,1]/m[i]
140 e = mE[i]/m[i]-fhalf * (u[i, 0]**2 + u[i, 1]**2)
141 p[i] = pfunc(e, rho[i], fluidlaw_data)
142 c[i] = cfunc(p[i], rho[i], fluidlaw_data)
144 def acoustic_impedance_tensor(self, n, outZ = None):
145 crho, cc = self.get_primal_variable(["rho", "c"])
146 if outZ is None :
147 batch_out = np.broadcast_shapes(c.shape, n.shape[:-1])
148 out_dtype = n.dtype
149 outZ = np.empty(batch_out + (DIM, DIM), dtype=out_dtype)
150 # if len(n.shape[:-1]) != len(crho.shape):
151 # crho = crho[:, np.newaxis]
152 # cc = cc[:, np.newaxis]
153 #_acoustic_impedance_tensor_2d_guvec(crho, cc, n, outZ)
154 _acoustic_impedance_tensor_2d_njit(crho, cc, n, outZ)
155 return outZ
157def GELEES_bc_set_sliding_wall(
158 vids: IndexArray,
159 NxN: FloatArray,
160 ImNxN: FloatArray,
161 vA: FloatArray,
162 vb: FloatArray
163) -> None:
164 r"""
165 Applies a sliding wall (slip) boundary condition to the global nodal
166 linear system.
168 This function modifies the assembled $A_r$ matrices and $b_r$ vectors
169 in-place for the nodes located on a sliding wall. It ensures that the
170 velocity vector resulting from the linear system solver will have zero
171 component normal to the wall ($u \cdot n = 0$).
173 The algebraic modification projects the system onto the tangent plane:
174 $A_{new} = (I - n \otimes n) A_r (I - n \otimes n) + \mathrm{Tr}(A_r)(n \otimes n)$
175 $b_{new} = (I - n \otimes n) b_r$
177 Adding the trace back to the normal component ensures the modified
178 matrix remains strictly invertible and well-conditioned.
180 Parameters
181 ----------
182 vids : npt.NDArray[np.uint64]
183 1D array of length `(n_boundary_nodes)`. The global indices of the
184 vertices located on the sliding wall.
185 NxN : FloatArray
186 3D array of shape `(n_boundary_nodes, 2, 2)`. The dyadic tensor
187 product of the normal vectors ($n \otimes n$) at each boundary node.
188 ImNxN : FloatArray
189 3D array of shape `(n_boundary_nodes, 2, 2)`. The tangent projection
190 matrix ($I - n \otimes n$) at each boundary node.
191 vA : FloatArray
192 3D array of shape `(nverts, 2, 2)`. The global assembled $A_r$
193 matrices for all nodes in the mesh. Modified in-place.
194 vb : FloatArray
195 2D array of shape `(nverts, 2)`. The global assembled $b_r$ vectors
196 for all nodes in the mesh. Modified in-place.
198 Returns
199 -------
200 None
201 The global `vA` and `vb` arrays are modified in-place to reflect
202 the boundary constraints.
203 """
204 boundary_vA = vA[vids]
205 # Extract the trace (A00 + A11) for the boundary nodes
206 trace_boundary_vA = boundary_vA[..., 0,0] + boundary_vA[..., 1,1]
207 # Apply the projection and add the penalized normal component
208 vA[vids] = np.matmul(np.matmul(ImNxN, boundary_vA), ImNxN) \
209 + trace_boundary_vA[...,np.newaxis, np.newaxis]*NxN
210 # Project the RHS vector onto the tangent plane
211 vb[vids] = np.matvec(ImNxN, vb[vids] )
213def GELEES_bc_set_external_pressure(
214 vids: IndexArray,
215 ps: FloatArray,
216 mesh: GELEES_MeshStruct,
217 out_vA: FloatArray, # pylint: disable=w0613
218 out_vb: FloatArray
219) -> None:
220 r"""Applies external pressure directly to the RHS vectors $b_{jr}$ of boundary nodes."""
221 for vid, p in zip(vids, ps) :
222 s = mesh.v2t_index[vid]
223 e = mesh.v2t_index[vid+1]
224 for cell, cv in zip(mesh.v2t_tris[s:e], mesh.v2t_v[s:e]):
225 out_vb[vid] -= p*mesh.C[cell, cv]
227FType_v = None
228FType_vc = None
230if USE_NUMBA:
232 @numba.njit(parallel=True, cache =True)
233 def rho_u_p_c_numba(U, S, pfunc, cfunc, fluidlaw_data, out_W):
234 N = U.shape[0]
235 for i in numba.prange(N): # pylint: disable=E1133
236 m_i, mu_i, me_i = get_conservative2d(i, U)
237 S_i = S[i]
238 rho_i = m_i/S_i
239 u_i = vec2_scal(1./m_i, mu_i)
240 e_i = me_i / m_i - npfloat(0.5) * vec2_square(u_i)
241 p_i = pfunc(e_i, rho_i, fluidlaw_data)
242 c_i = cfunc(p_i, rho_i, fluidlaw_data)
243 set_primal2d(rho_i, u_i, p_i, c_i, i, out_W)
245 # @numba.njit(cache=True, fastmath = FASTMATH)
246 # def GELEES_tri_dt(Sj, nCj0, nCj1, nCj2, cj):
247 # ''' return th critical time step in a triangle'''
248 # sum_normC = nCj0 + nCj1 + nCj2
249 # dt = Sj / sum_normC/ cj
250 # return dt
252 # @numba.njit(parallel=True, cache=True, fastmath = FASTMATH)
253 # def critical_timestep_numba(S, nC, c):
254 # min_dt = np.inf # numba does autoamatique reduction on this variable
255 # for j in numba.prange(len(S)): # pylint: disable=E1133
256 # Sj = S[j]
257 # nCj0 = nC[j,0]
258 # nCj1 = nC[j,1]
259 # nCj2 = nC[j,2]
260 # cj = c[j]
261 # dt = GELEES_tri_dt(Sj, nCj0, nCj1, nCj2, cj)
262 # # # Numba comprend automatiquement qu'il s'agit d'une réduction parallèle
263 # min_dt = min(min_dt, dt)
264 # return min_dt
266 @numba.njit(parallel=False, cache=True, fastmath = FASTMATH)
267 def GELEES_getAbj(rhoj:npfloat, uj:Vec2, pj:npfloat, cj:npfloat, Cjr:Vec2)->Tuple[Sy22,Vec2]:
268 ''' compute Ajr and bjr which ar part of the node k of cellj
269 linear systeme to solve at node jr '''
270 # Aj = rhoj*cj/|Cjr|*cjr X cjr
271 # bj = pj*Cjr + Aj.uj
272 scalj = rhoj * cj/vec2_norm(Cjr)
273 Aj = sy22_r(scalj, Cjr)
274 Auj = sy22_mv(Aj, uj)
275 bj = vec2_add( vec2_scal(pj, Cjr), Auj)
276 return Aj, bj
278 @numba.njit(cache=True, fastmath = FASTMATH)
279 def bc_set_sliding_wall_numba(n, A, b):
280 An = sy22_mv(A, n)
281 scal = vec2_dot(n, An) + sy22_trace(A)
282 A = sy22_add(A, sy22_r(scal, n))
283 A = sy22_add(A, sy22_r2(npfloat(-fone), n, An))
284 b = vec2_add(b, vec2_scal(-vec2_dot(n, b), n))
285 return A, b
287 def as_numba_pointer(func, signature):
288 """
289 Encapsule une fonction Python dans une liste typée Numba.
290 C'est un hack nécessaire pour forcer Numba à utiliser un pointeur
291 générique (First-Class Function) et éviter de casser le cache JIT.
292 """
293 # type_list_of_func = nbt.ListType(signature)
294 # # On instancie la liste normalement avec l'argument caché 'lsttype'
295 # ptr_box = nbt.List(dtype= type_list_of_func)
296 # ptr_box.append(func)
299 creator = getattr(nbty.List, "empty_list")
301 # 2. On instancie le conteneur physique en lui passant la signature
302 ptr_box = creator(signature)
304 # 3. Maintenant, .append() fonctionnera parfaitement !
305 ptr_box.append(func)
306 return ptr_box
309 data_1d_t = nbfloat[:]
310 # Signature pour les Nœuds : (Ar, br) = f(Ar, br, data_1d)
311 signature_vbc = nbct.Tuple((sy22_t, vec2_t))(sy22_t, vec2_t, data_1d_t)
313 # Signature pour les Coins : (Aj, bj) = f(Aj, bj, Cjr, data_1d)
314 signature_vcbc = nbct.Tuple((sy22_t, vec2_t))(sy22_t, vec2_t, vec2_t, data_1d_t)
316 # On crée les "Types" de Fonction Abstraite
317 FType_v = nbct.FunctionType(signature_vbc)
318 FType_vc = nbct.FunctionType(signature_vcbc)
320 @numba.njit(signature_vbc, cache=True, fastmath=FASTMATH)
321 def bc_vnone_cache(Ar, br, data): # pylint: disable=W0613
322 return Ar, br
324 @numba.njit(signature_vbc, cache=True, fastmath=FASTMATH)
325 def bc_sliding_wall_cache(Ar, br, data): # pylint: disable=w0613
326 n = (data[0], data[1])
327 return bc_set_sliding_wall_numba(n, Ar, br)
329 @numba.njit(signature_vcbc, cache=True, fastmath=FASTMATH)
330 def bc_vcnone_cache(Aj, bj, Cjr, data): # pylint: disable=w0613
331 return Aj, bj
333 @numba.njit(signature_vcbc, cache=True, fastmath=FASTMATH)
334 def bc_pressure_cache(Aj, bj, Cjr, data):
335 p = data[0]
336 bj_new = vec2_sub(bj, vec2_scal(p, Cjr))
337 return Aj, bj_new
339 @numba.njit(inline='always')
340 def apply_bc_v(ptr_box, Ar, br, data):
341 """ Cache la misère du [0] pour les noeuds """
342 f = ptr_box[0]
343 return f(Ar, br, data)
345 @numba.njit(inline='always')
346 def apply_bc_vc(ptr_box, Aj, bj, Cjr, data):
347 """ Cache la misère du [0] pour les coins """
348 f = ptr_box[0]
349 return f(Aj, bj, Cjr, data)
351 def GELEES_bulk_vertices_loop_numba(vids, mesh, W, dt, xy, vu):
352 return _GELEES_bulk_vertices_loop_numba(vids, mesh, W, dt, xy, vu)
354 @numba.njit(parallel=True, cache=True, fastmath = FASTMATH)
355 def _GELEES_bulk_vertices_loop_numba(vids, mesh, W, dt, out_xy, out_vu):
356 # for each vertex in vids, construct the 2x2 linear system
357 # to find the vertex velocity and update the node position
358 for i in numba.prange(len(vids)): # pylint: disable=E1133
359 vid = vids[i]
360 # setup riman solver at node i
361 s = mesh.v2t_index[vid]
362 e = mesh.v2t_index[vid+1]
363 Ar = sy22_zero
364 br = vec2_zero
365 # Cell k connectd to node vid
366 for k in range(s, e):
367 tj = mesh.v2t_tris[k]
368 vr = mesh.v2t_v[k]
369 Cjr = (mesh.C[tj, vr, 0], mesh.C[tj, vr, 1])
370 rhoj, uj, pj, cj = get_primal2d(tj, W)
371 Aj, bj = GELEES_getAbj(rhoj, uj, pj, cj, Cjr)
372 Ar = sy22_add(Ar, Aj)
373 br = vec2_add(br, bj)
374 # solve rieman problem (find correct velocity)
375 vur = sy22_smv(Ar, br)
376 set_vec2(vid, out_vu, vur)
377 # update node position
378 xyr = vec2_add(get_vec2(vid, out_xy), vec2_scal(dt, vur))
379 set_vec2(vid, out_xy, xyr)
381 def GELEES_specialized_vertices_loop(vids, mesh, W, dt, box_f_v, box_f_vc, data_v, data_vc,\
382 out_xy, out_vu):
383 _GELEES_specialized_vertices_loop(vids, mesh, W, dt, box_f_v, box_f_vc, data_v, data_vc,\
384 out_xy, out_vu)
386 @numba.njit(parallel=True, cache=True, fastmath = FASTMATH)
387 def _GELEES_specialized_vertices_loop(vids, mesh, W, dt, box_f_v, box_f_vc, data_v, data_vc,\
388 out_xy, out_vu):
389 # for each vertex in vids, construct the 2x2 linear system
390 # to find the vertex velocity and update the node position
391 for i in numba.prange(len(vids)): # pylint: disable=E1133
392 vid = vids[i]
393 s = mesh.v2t_index[vid]
394 e = mesh.v2t_index[vid+1]
395 Ar = sy22_zero # matrice 22 symetrique stockée ainsi (A00, A01, A11)
396 br = vec2_zero
397 node_data_v = data_v[i]
398 node_data_vc = data_vc[i]
399 # Cell k connectd to node vid
400 for k in range(s, e):
401 tj = mesh.v2t_tris[k]
402 vr = mesh.v2t_v[k]
403 Cjr = (mesh.C[tj, vr, 0], mesh.C[tj, vr, 1])
404 rhoj, uj, pj, cj = get_primal2d(tj, W)
405 Aj, bj = GELEES_getAbj(rhoj, uj, pj, cj, Cjr)
406 Aj, bj = apply_bc_vc(box_f_vc, Aj, bj, Cjr, node_data_vc)
407 Ar = sy22_add(Ar, Aj)
408 br = vec2_add(br, bj)
409 Ar, br = apply_bc_v(box_f_v, Ar, br, node_data_v)
410 vur = sy22_smv(Ar, br)
411 set_vec2(vid, out_vu, vur)
412 xyr = vec2_add(get_vec2(vid, out_xy), vec2_scal(dt, vur))
413 set_vec2(vid, out_xy, xyr)
416 @numba.njit(parallel=False, cache=True, fastmath = FASTMATH)
417 def cell_vertex_pressure(cp, crho, cc, cu, vu, vN):
418 return cp - crho*cc*vec2_dot(vec2_sub(vu,cu), vN)
420 def GELEES_cell_loop(dt, msh, vu, U, W, pfunc, cfunc, fluidlaw_data):
421 return _GELEES_cell_loop(dt, msh, vu, U, W, pfunc, cfunc, fluidlaw_data)
423 @numba.njit(parallel=True, cache=True, fastmath=FASTMATH)
424 def _GELEES_cell_loop(dt: npfloat, mesh, vu: FloatArray, U: FloatArray, W: FloatArray,\
425 pfunc, cfunc, fluidlaw_data) -> float:
426 r"""
427 Ultra-optimized cell loop: Fully unrolled (no inner loops),
428 manual prefetching to registers, minimized divisions.
429 """
430 min_dtc = np.inf
431 # Extraction des pointeurs bruts au niveau global
432 t2v, nC, C, S, xy = mesh.t2v, mesh.nC, mesh.C, mesh.cell_meas, mesh.xy
434 for j in numba.prange(mesh.ncells): # pylint: disable=E1133
435 # =========================================================
436 # 1. PREFETCHING (Lecture RAM -> Registres L1)
437 # =========================================================
438 n0, n1, n2 = t2v[j, 0], t2v[j, 1], t2v[j, 2]
439 vur0, vur1, vur2 = get_vec2(n0, vu), get_vec2(n1, vu), get_vec2(n2, vu)
440 Cj0, Cj1, Cj2 = get_tri_vec2(j, C)
442 # Thermo courante
443 mj, muj, mej = get_conservative2d(j, U)
444 rhoj, uj, pj, cj = get_primal2d(j, W)
446 # =========================================================
447 # 2. CALCUL DES FLUX (Désenroulé, vectorisé SIMD)
448 # =========================================================
449 # Inverses des normes (3 divisions factorisées)
450 nC0, nC1, nC2 = get_tri_scal(j, nC)
451 Nj0, Nj1, Nj2 = vec2_scal(fone/nC0, Cj0),\
452 vec2_scal(fone/nC1, Cj1), \
453 vec2_scal(fone/nC2, Cj2)
454 # Pressions nodales (Résolveur de Riemann acoustique local)
455 pjr0 = cell_vertex_pressure(pj, rhoj, cj, uj, vur0, Nj0)
456 pjr1 = cell_vertex_pressure(pj, rhoj, cj, uj, vur1, Nj1)
457 pjr2 = cell_vertex_pressure(pj, rhoj, cj, uj, vur2, Nj2)
458 # Quantité de mouvement : q_j = sum(pjr * Cjr)
459 qfluxj = vec2_add3( vec2_scal(pjr0, Cj0), vec2_scal(pjr1, Cj1), vec2_scal(pjr2, Cj2))
460 # Énergie : e_j = sum(pjr * dot(vur, Cjr))
461 efluxj = pjr0*vec2_dot(vur0, Cj0) + pjr1*vec2_dot(vur1, Cj1) + pjr2*vec2_dot(vur2, Cj2)
462 # =========================================================
463 # 3. MISE A JOUR CONSERVATIVE
464 # =========================================================
465 muj_new = vec2_sub(muj, vec2_scal(dt, qfluxj) )
466 mej_new = mej - dt * efluxj
467 # =========================================================
468 # 4. NOUVELLE GEOMETRIE (On réutilise n0, n1, n2 !)
469 # =========================================================
470 xy0, xy1, xy2 = get_vec2(n0, xy), get_vec2(n1, xy), get_vec2(n2, xy)
471 Sj_new, Cj0_new, Cj1_new, Cj2_new = Cjr_triangle(xy0, xy1, xy2)
472 nCj0_new = vec2_norm(Cj0_new)
473 nCj1_new = vec2_norm(Cj1_new)
474 nCj2_new = vec2_norm(Cj2_new)
475 # =========================================================
476 # 5. NOUVELLE THERMODYNAMIQUE
477 # =========================================================
478 inv_Sj_new = 1.0/Sj_new
479 rhoj_new = mj * inv_Sj_new
480 inv_mj = 1.0 / mj
481 uj_new = vec2_scal(inv_mj, muj_new)
482 ej_new = mej_new * inv_mj - fhalf * vec2_square(uj_new)
483 pj_new = pfunc(ej_new, rhoj_new, fluidlaw_data)
484 cj_new = cfunc(pj_new, rhoj_new, fluidlaw_data)
485 # =========================================================
486 # 6. ECRITURE MEMOIRE ET PAS DE TEMPS
487 # =========================================================
488 set_conservative2d(mj, muj_new, mej_new, j, U)
489 set_primal2d(rhoj_new, uj_new, pj_new, cj_new, j, W)
490 S[j] = Sj_new
491 set_tri_scal(j, nC, nCj0_new, nCj1_new, nCj2_new)
492 set_tri_vec2(j, C, Cj0_new, Cj1_new, Cj2_new)
493 # Pas de temps critic
494 dtc = Sj_new / ((nCj0_new + nCj1_new + nCj2_new) * cj_new)
495 min_dtc = min(min_dtc, dtc)
496 return min_dtc
498 @numba.njit(inline='always', fastmath = FASTMATH)
499 def get_vec2(index, array):
500 return (array[index, 0], array[index, 1])
502 @numba.njit(cache=True, fastmath = FASTMATH)
503 def get_tri_vec2(index, array):
504 return (array[index, 0, 0], array[index,0, 1]), \
505 (array[index, 1, 0], array[index,1, 1]), \
506 (array[index, 2, 0], array[index,2, 1]),
508 @numba.njit(inline='always', fastmath = FASTMATH)
509 def get_tri_scal(index, array):
510 return array[index, 0], array[index, 1], array[index, 2]
512 @numba.njit(inline='always', fastmath = FASTMATH)
513 def set_vec2(index, array, v):
514 array[index, 0], array[index, 1] = v[0], v[1]
516 @numba.njit(inline='always', fastmath = FASTMATH)
517 def set_tri_scal(index, array, s0, s1, s2):
518 array[index, 0], array[index, 1], array[index, 2] = s0, s1, s2
520 @numba.njit(cache= True, fastmath = FASTMATH)
521 def set_tri_vec2(index, array, v0, v1, v2):
522 array[index, 0, 0], array[index, 0, 1] = v0[0], v0[1]
523 array[index, 1, 0], array[index, 1, 1] = v1[0], v1[1]
524 array[index, 2, 0], array[index, 2, 1] = v2[0], v2[1]
526 @numba.njit(cache=True, fastmath = FASTMATH)
527 def get_trixy(elemid, mesh):
528 tri = mesh.t2v[elemid]
529 n0, n1, n2 = tri[0], tri[1], tri[2]
530 xy = mesh.xy
531 xy0 = get_vec2(n0, xy)
532 xy1 = get_vec2(n1, xy)
533 xy2 = get_vec2(n2, xy)
534 return xy0, xy1, xy2
536 @numba.njit(cache=True, fastmath = FASTMATH)
537 def get_conservative2d(index, U):
538 m, mu, me = U[index, I_M], (U[index, I_MU_START], U[index, I_MU_START+1]), U[index, I_ME]
539 return m, mu, me
541 @numba.njit(cache=True, fastmath = FASTMATH)
542 def set_conservative2d(m, mu, me, index, U):
543 U[index, I_M], U[index, I_MU_START], U[index, I_MU_START+1], U[index, I_ME] = m, mu[0], mu[1], me
545 @numba.njit(cache=True, fastmath = FASTMATH)
546 def get_primal2d( index, W):
547 rho, u, p, c = W[index,I_RHO], (W[index,I_U_START], W[index,I_U_START+1]), W[index,I_P], W[index,I_C]
548 return rho, u, p, c
550 @numba.njit(cache=True, fastmath = FASTMATH)
551 def set_primal2d(rho, u, p, c, index, W):
552 W[index, I_RHO], W[index, I_U_START], W[index, I_U_START+1], W[index, I_P], W[index, I_C] = rho, u[0], u[1], p, c
557class GlaceExplicitLagrangian_EulerEquationSolver(GlaceExplicitLagrangianBaseSolver):# pylint: disable=C0103
558 """
559 A 2D Explicit Lagrangian Finite Volume Solver for the Euler Equations
560 using the GLACE (Godunov-type Lagrangian Scheme for Compressible Euler) scheme.
562 This class manages the mesh topology, physical state variables, boundary
563 conditions, and time integration (both in pure NumPy and Numba JIT modes).
564 """
565 def __init__(
566 self,
567 xy: FloatArray,
568 tris: IndexArray,
569 boundary_conditions: Dict[str, Any],
570 use_numba: bool = True
571 ) -> None:
572 """
573 Initializes the GLACE solver.
575 Parameters
576 ----------
577 xy : FloatArray
578 2D array of shape `(nverts, 2)` containing the initial nodal coordinates.
579 tris : IndexArray
580 2D array of shape `(ncells, 3)` containing the cell-to-vertex connectivity.
581 boundary_conditions : Dict[str, Any]
582 Dictionary containing the boundary condition definitions (types, nodes,
583 normals, and prescribed pressures).
584 use_numba : bool, optional
585 If True, the solver will attempt to use the optimized Numba JIT-compiled
586 routines. Defaults to True.
587 """
588 self.use_numba = False
589 if use_numba :
590 if USE_NUMBA:
591 self.use_numba = True
592 else :
593 print("Warning : Trying to use numba will the __FILE__ USE_NUMBA is set to false")
594 self.gaz = PerfectGaz()
596 mesh = setup_MeshStruct(xy, tris)
597 nverts = mesh.nverts
598 ncells = mesh.ncells
599 conservation_law = ConservationLawManager_Euler(self.gaz, ncells)
600 super().__init__(conservation_law, mesh)
601 self.t = 0.
602 self.vu = np.zeros((nverts, DIM), dtype=npfloat)
603 # conservative vars (m, mux, muy, me) per cell
604 # primary vars (rho, ux, uy, p, c) per cell
605 self._init_bc(boundary_conditions)
607 if self.use_numba :
608 #computeCjr_triangles_numba(self.mesh)
609 self.numba_pressure_c_func, self.numba_soundspeed_c_func, self.fluidlaw_data = self.gaz.get_numba_c_funcs()
610 self.step = self.step_numba
611 else:
612 #Cjr_triangle_np(self.mesh)
613 self.step = self.step_np
615 def _init_bc(self, boundary_conditions: Dict[str, Any]) -> None:
616 # 1. ACCUMULATION : On récolte les infos physiques par noeud
617 # Chaque noeud a un dictionnaire avec ses normales et ses pressions
618 node_data = defaultdict(lambda: {'normals': [], 'pressures': []})
620 for bc_data in boundary_conditions.values():
621 bc_type = bc_data["type"]
622 vids = bc_data["vid"]
624 if bc_type == "sliding_wall":
625 n = np.array(bc_data["n"], dtype=npfloat)
626 for vid in vids:
627 node_data[vid]['normals'].append(n)
628 elif bc_type == "external_pressure":
629 p = float(bc_data["pressure"])
630 for vid in vids:
631 node_data[vid]['pressures'].append(p)
633 # 2. CLASSIFICATION : Préparation des listes par catégorie
634 categories = {
635 "sliding_wall": {"nodes": [], "N": []},
636 "external_pressure": {"nodes": [], "pressures": []},
637 "pressure_and_wall": {"nodes": [], "pressures": [], "N": []},
638 "fixed_point": {"nodes": []}
639 }
641 for vid, data in node_data.items():
642 normals = data['normals']
643 pressures = data['pressures']
644 n_walls = len(normals)
645 n_press = len(pressures)
646 if n_walls == 1 and n_press == 0:
647 categories["sliding_wall"]["nodes"].append(vid)
648 categories["sliding_wall"]["N"].append(normals[0])
649 elif n_walls == 0 and n_press >= 1:
650 categories["external_pressure"]["nodes"].append(vid)
651 categories["external_pressure"]["pressures"].append(np.mean(pressures))
652 elif n_walls == 1 and n_press >= 1:
653 categories["pressure_and_wall"]["nodes"].append(vid)
654 categories["pressure_and_wall"]["N"].append(normals[0])
655 categories["pressure_and_wall"]["pressures"].append(np.mean(pressures))
656 elif n_walls >= 2:
657 # 2 murs ou plus qui se croisent = point fixe (coin)
658 categories["fixed_point"]["nodes"].append(vid)
659 else:
660 print(f"Warning: BC non gérée pour le noeud {vid}")
662 # 3. VECTORISATION : Création du dictionnaire final avec les arrays NumPy
663 self.bc = {}
664 for cat_name, cat_data in categories.items():
665 bc_dict ={} #{"nodes": np.array(cat_data["nodes"], dtype=IndexType)}
666 bc_dict["nodes"] = np.array(cat_data["nodes"], dtype=IndexType)
667 if "N" in cat_data:
668 N_arr = np.array(cat_data["N"], dtype=npfloat).reshape(-1, 2)
669 bc_dict["N"] = N_arr
670 # Pour le solveur NumPy pur, on vectorise la création des tenseurs (NxN et I-NxN)
671 if not self.use_numba:
672 if len(N_arr) > 0:
673 NxN = N_arr[:, :, np.newaxis] * N_arr[:, np.newaxis, :]
674 bc_dict["NxN"] = NxN
675 bc_dict["ImNxN"] = np.eye(2) - NxN
676 else:
677 # Ensures the keys exist even if there are no nodes
678 bc_dict["NxN"] = np.empty((0, 2, 2), dtype=npfloat)
679 bc_dict["ImNxN"] = np.empty((0, 2, 2), dtype=npfloat)
681 if "pressures" in cat_data:
682 bc_dict["pressures"] = np.array(cat_data["pressures"], dtype=npfloat)
683 self.bc[cat_name] = bc_dict
684 # Separating Bulk from boundary vertices
685 boundary_vids = np.concatenate([bc["nodes"] for bc in self.bc.values()])
686 self.bulk_vids = np.setdiff1d(np.arange(self.mesh.nverts), boundary_vids)
687 # BINDING NUMBA : attaching function pointer and data for each bc type.
688 if self.use_numba:
689 for bc_type, bc_dict in self.bc.items():
690 n_nodes = len(bc_dict["nodes"])
691 if bc_type == "sliding_wall":
692 bc_dict["vdata_cache"] = bc_dict["N"]
693 bc_dict["vfun_cache"] = as_numba_pointer(bc_sliding_wall_cache, FType_v)
694 bc_dict["vcdata_cache"] = np.empty((n_nodes, 0), dtype=npfloat)
695 bc_dict["vcfun_cache"] = as_numba_pointer(bc_vcnone_cache, FType_vc)
696 elif bc_type == "external_pressure":
697 bc_dict["vdata_cache"] = np.empty((n_nodes, 0), dtype=npfloat)
698 bc_dict["vfun_cache"] = as_numba_pointer(bc_vnone_cache, FType_v)
699 bc_dict["vcdata_cache"] = bc_dict.get("pressures", np.empty(0, dtype=npfloat)).reshape(-1, 1)
700 bc_dict["vcfun_cache"] = as_numba_pointer(bc_pressure_cache, FType_vc)
701 elif bc_type == "pressure_and_wall":
702 bc_dict["vdata_cache"] = bc_dict["N"]
703 bc_dict["vfun_cache"] = as_numba_pointer(bc_sliding_wall_cache, FType_v)
704 bc_dict["vcdata_cache"] = bc_dict.get("pressures", np.empty(0, dtype=npfloat)).reshape(-1, 1)
705 bc_dict["vcfun_cache"] = as_numba_pointer(bc_pressure_cache, FType_vc)
706 # Nothing done for "fixed_point" : we just have to make sure those points are never moved.
708 def _set_boundary_conditions(self)-> None:
709 # external_pressure
710 bc_data = self.bc["external_pressure"]
711 vids = bc_data["nodes"]
712 ps = bc_data["pressures"]
713 GELEES_bc_set_external_pressure(vids, ps, self.mesh, self.vA, self.vb)
714 # sliding wall :
715 bc_data = self.bc["sliding_wall"]
716 vids = bc_data["nodes"]
717 NxN = bc_data["NxN"]
718 ImNxN = bc_data["ImNxN"]
719 GELEES_bc_set_sliding_wall(vids, NxN, ImNxN, self.vA, self.vb)
720 # pressure_and_wall
721 bc_data = self.bc["pressure_and_wall"]
722 vids = bc_data["nodes"]
723 ps = bc_data["pressures"]
724 NxN = bc_data["NxN"]
725 ImNxN = bc_data["ImNxN"]
726 GELEES_bc_set_external_pressure(vids, ps, self.mesh, self.vA, self.vb)
727 GELEES_bc_set_sliding_wall(vids, NxN, ImNxN, self.vA, self.vb)
728 # Fixed (2 sliding walls)
729 bc_data = self.bc["fixed_point"]
730 vids = bc_data["nodes"]
731 self.vb[vids] = 0.
733 def _riemann_solver_setup(self)-> None:
734 cvn: FloatArray = self.mesh.N
735 cvnC: FloatArray = self.mesh.nC
736 cvC: FloatArray = self.mesh.C
737 cvA: FloatArray = self.cvA
738 cvb: FloatArray = self.cvb
739 cu, cp = self.conservation_law.get_primal_variable(["u", "p"])
740 self.conservation_law.acoustic_impedance_tensor(cvn, outZ = self.cvA)
741 if USE_NUMBA :
742 GlaceExplicitLagrangian_EulerEquationSolver._riemann_solver_setup_njit(cu, cp, cvA, cvC, cvnC, cvA, cvb)
743 else:
744 cvA *= cvnC[..., np.newaxis, np.newaxis]
745 cp_ = np.expand_dims(cp, axis = (1,2))
746 # b_jr = p * Cjr + A_jr * u
747 cvb[:] = cvC*cp_ \
748 + self.cvA[...,0]*cu[:, np.newaxis, np.newaxis, 0] \
749 + self.cvA[...,1]*cu[:, np.newaxis, np.newaxis, 1]
751 if USE_NUMBA :
752 @staticmethod
753 @numba.njit(parallel=True, cache = True, fastmath = FASTMATH)
754 def _riemann_solver_setup_njit(cu, cp, cvZ, cvC, cvnC, out_cvA, out_cvb):
755 N =cu.shape[0]
756 for j in numba.prange(N): # type ignore # pylint: disable=E1133
757 pj = cp[j]
758 uj0 = cu[j,0]
759 uj1 = cu[j,1]
760 vAj = out_cvA[j]
761 vZj = cvZ[j]
762 vnCj = cvnC[j]
763 vbj = out_cvb[j]
764 vCj = cvC[j]
765 for r in range(3):
766 Cjr = vCj[r]
767 nCjr = vnCj[r]
768 Zjr = vZj[r]
769 Ajr = vAj[r]
770 Ajr[0, 0] = nCjr*Zjr[0 ,0]
771 Ajr[0, 1] = nCjr*Zjr[0 ,1]
772 Ajr[1, 0] = nCjr*Zjr[1 ,0]
773 Ajr[1, 1] = nCjr*Zjr[1 ,1]
774 bjr = vbj[r]
775 bjr[0] = pj*Cjr[0] + Ajr[0, 0]*uj0 + Ajr[0, 1]*uj1
776 bjr[1] = pj*Cjr[1] + Ajr[1, 0]*uj0 + Ajr[1, 1]*uj1
778 def set_initial_state(
779 self,
780 density: Callable[[FloatArray], FloatArray],
781 velocity: Callable[[FloatArray], FloatArray],
782 pressure: Callable[[FloatArray], FloatArray]
783 ) -> None:
784 """
785 Sets the initial macroscopic state of the fluid for all cells.
787 Parameters
788 ----------
789 density : Callable[[FloatArray], FloatArray]
790 A function taking cell center coordinates (N, 2) and returning
791 an array of initial cu, cp = self.conservation_law.get_primal_variable(["u", "p"])
792 cp_ = np.expand_dims(cp, axis = (1,2))
793 # b_jr = p * Cjr + A_jr * u
794 self.cvb[:] = cvC*cp_ \
795 + self.cvA[...,0]*cu[:, np.newaxis, np.newaxis, 0] \
796 + self.cvA[...,1]*cu[:, np.newaxis, np.newaxis, 1]
797densities (N,).
798 velocity : Callable[[FloatArray], FloatArray]
799 A function taking cell center coordinates (N, 2) and returning
800 an array of initial velocity vectors (N, 2).
801 pressure : Callable[[FloatArray], FloatArray]
802 A function taking cell center coordinates (N, 2) and returning
803 an array of initial pressures (N,).
804 """
805 xytris = self.mesh.xy[self.mesh.t2v]
806 xycog = np.sum(xytris, axis =1)/3.
807 densitycog = np.array(density(xycog), dtype = npfloat)
808 velocitycog = np.array(velocity(xycog), dtype = npfloat)
809 pressurecog = np.array(pressure(xycog), dtype = npfloat)
810 cell_measure = self.mesh.cell_meas
811 m = densitycog*cell_measure
812 int_spec_enrg = self.gaz.internal_specific_energy(pressurecog,densitycog)
813 kin_spec_enrg = fhalf*(velocitycog[:, 0]**2 + velocitycog[:, 1]**2)
814 getU = self.conservation_law.get_conservative_variable
815 getU("m")[:] = m
816 getU("mu")[:] = velocitycog*m[..., np.newaxis]
817 getU("mE")[:] = (int_spec_enrg+kin_spec_enrg)*m
818 self.conservation_law.conservative_to_primal(cell_measure)
820 def critical_timestep(self) -> float:
821 """
822 Computes the maximum allowable critical time step for the entire mesh
823 based on the CFL condition (acoustic wave speeds).
825 Returns
826 -------
827 dtmax : float
828 The critical time step (dt).
829 """
830 V = self.mesh.cell_meas
831 cc = self.get_primal_variable("c") # per cell sound speed
832 nC = self.mesh.nC
833 return type(self).critical_timestep_pnjit(V, nC, cc)
835 # if True : #self.use_numba:
836 # nC = self.mesh.nC
837 # return critical_timestep_numba(V, nC, cc)
838 # else:
839 # C = self.mesh.C
840 # cvnormC = np.sqrt( np.square(C[...,0]) + np.square(C[...,1]) )
841 # #sum_cvnormC =np.sum(cvnormC, axis =1)
842 # sum_cvnormC = cvnormC[:,0] + cvnormC[:,1] + cvnormC[:,2]
843 # cL = V/sum_cvnormC
844 # dtmax = np.min(cL/cc)
845 # return dtmax
847 def step(self, dt: float) -> float: # pylint: disable=E0202:method-hidden
848 """
849 Advances the simulation by one time step `dt`.
851 Notes
852 -----
853 This method is dynamically routed at runtime. Depending on the `use_numba`
854 flag provided during the solver's initialization, calling `step()` will
855 execute either:
856 - `step_numba()` : The heavily optimized JIT-compiled solver.
857 - `step_np()` : The pure vectorized NumPy solver.
859 Parameters
860 ----------
861 dt : float
862 The duration of the time step.
864 Returns
865 -------
866 dtnew : float
867 The new computed critical time step for the next iteration.
868 """
869 # This line will never be executed thanks to dynamic affctation to self.step in __init__
870 raise NotImplementedError("This method is dynamically assigned during __init__.")
872 def step_np(self, dt: float) -> float:
873 """
874 Advances the simulation by one time step `dt` using the pure NumPy solver.
876 Parameters
877 ----------
878 dt : float
879 The duration of the time step.
881 Returns
882 -------
883 dtcnew : float
884 The new computed critical time step for the next iteration.
885 """
886 #nc = self.ncells
887 vu, cvqflux = self._riemann_solver()
888 mesh = self.mesh
889 cvu = self.cvu
890 #graph.extract_txy(vu, self.mesh.t2v, cvu)
891 #cvu = vu[self.mesh.t2v]
895 GlaceExplicitLagrangian_EulerEquationSolver.update_conservative(cvqflux, cvu, dt, self.conservation_law)
897 #moving node
898 mesh.xy[:] += dt*vu
899 # update primal variable (on moved cells)
900 update_MeshStruct(mesh.xy, mesh)
901 self.conservation_law.conservative_to_primal(mesh.cell_meas)
902 #update time
903 self.t += dt
904 #compute critical time step
905 dtcnew = self.critical_timestep()
906 return dtcnew
908 @staticmethod
909 def update_conservative(cvqflux, cvu, dt, cons_law_man):
910 mu, mE = cons_law_man.get_conservative_variable(["mu", "mE"])
911 if USE_NUMBA:
912 GlaceExplicitLagrangian_EulerEquationSolver.update_conservative_njit(cvqflux, cvu, dt, mu, mE)
913 else:
914 mu[:, 0] -= dt*(cvqflux[:,0,0] + cvqflux[:,1,0] + cvqflux[:,2,0])
915 mu[:, 1] -= dt*(cvqflux[:,0,1] + cvqflux[:,1,1] + cvqflux[:,2,1])
916 mE -= dt*( cvqflux[:,0,0] * cvu[:,0,0] + cvqflux[:,0,1] * cvu[:,0,1] \
917 + cvqflux[:,1,0] * cvu[:,1,0] + cvqflux[:,1,1] * cvu[:,1,1] \
918 + cvqflux[:,2,0] * cvu[:,2,0] + cvqflux[:,2,1] * cvu[:,2,1])
919 if USE_NUMBA:
920 @staticmethod
921 @numba.njit(parallel = True, cache = True, fastmath = FASTMATH )
922 def update_conservative_njit(cvqflux, cvu, dt, mu, mE):
923 for j in numba.prange(cvqflux.shape[0]): # type ignore # pylint: disable=E1133
924 vqfluxj = cvqflux[j]
925 vuj = cvu[j]
926 mu[j, 0] -= dt*(vqfluxj[0,0] + vqfluxj[1,0] + vqfluxj[2,0])
927 mu[j, 1] -= dt*(vqfluxj[0,1] + vqfluxj[1,1] + vqfluxj[2,1])
928 mE[j] -= dt*( vqfluxj[0,0] * vuj[0,0] + vqfluxj[0,1] * vuj[0,1] \
929 + vqfluxj[1,0] * vuj[1,0] + vqfluxj[1,1] * vuj[1,1] \
930 + vqfluxj[2,0] * vuj[2,0] + vqfluxj[2,1] * vuj[2,1])
933 if USE_NUMBA :
934 def step_numba(self, dt: float) -> float:
935 """
936 Advances the simulation by one time step `dt` using the optimized
937 Numba JIT-compiled solver.
939 Parameters
940 ----------
941 dt : float
942 The duration of the time step.
944 Returns
945 -------
946 dtnew : float
947 The new computed critical time step for the next iteration.
948 """
949 mesh = self.mesh # the msh
950 W = self.conservation_law.W # cell primal variable
951 U = self.conservation_law.U # cell conservative variable
952 vu = self.vu # vertex velocity
953 # 1. bulk nodes update
954 bulk_vids = self.bulk_vids
955 #GELEES_vertices_loop(bulk_vids, mesh, W, dt, bc_vnone, bc_vcnone, mesh.xy, vu)
956 GELEES_bulk_vertices_loop_numba(bulk_vids, mesh, W, dt, mesh.xy, vu)
958 # 2. boundary node update
959 # 2.1 wall bc
961 for bc_type in ["sliding_wall", "external_pressure", "pressure_and_wall"]:
962 bc = self.bc[bc_type]
963 bc_vids = bc["nodes"]
964 if len(bc_vids) > 0 :
965 box_bc_vfun = bc["vfun_cache"]
966 box_bc_vcfun = bc["vcfun_cache"]
967 v_data = bc["vdata_cache"]
968 vc_data = bc["vcdata_cache"]
969 GELEES_specialized_vertices_loop(bc_vids, mesh, W, dt, box_bc_vfun, box_bc_vcfun, v_data, vc_data, mesh.xy, vu)
970 p_c_func = self.numba_pressure_c_func
971 sp_c_func = self.numba_soundspeed_c_func
972 fluidlaw_data = self.fluidlaw_data
973 dtnew = GELEES_cell_loop(dt, mesh, vu, U, W, p_c_func, sp_c_func, fluidlaw_data)
974 self.t += dt
975 return dtnew