#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# pylint: disable=C0103, W1401, C0302
r"""
Created on Thu Mar 7 11:17:17 2024
@author: nchevaug
2D Explicit Lagrangian Finite Volume Solver (GLACE Scheme).
Philosophy & Architecture:
This module implements the Godunov-type Lagrangian Scheme for Compressible
Euler equations (GLACE). It computes the evolution of a fluid on an unstructured
2D triangular mesh moving with the fluid velocity (Lagrangian framework).
The solver features a dual-execution architecture:
1. Pure NumPy (Fallback): Fully vectorized array operations.
2. Numba JIT (HPC Mode): Aggressively optimized loop-level parallelism
(`prange`), avoiding heap allocations and leveraging hardware registers
for the nodal Riemann solver.
"""
from collections import defaultdict
from typing import Any, Callable, Dict, Tuple, NamedTuple, Sequence, List, Union, assert_never, Optional, overload
from numba import guvectorize, float32, float64
from numba.np.ufunc import parallel
import numpy as np
from pyxfem.core.config import (
FASTMATH,
FloatArray,
IndexArray,
IndexType,
USE_NUMBA,
fhalf,
fone,
npfloat,
)
from pyxfem.fluidmechanics.exact_riemann_solver import PerfectGaz
from pyxfem.core.finite_volume_2d import (
ConservativeLawManager,
GELEES_MeshStruct,
setup_MeshStruct,
update_MeshStruct,
GlaceExplicitLagrangianBaseSolver,
Cjr_triangle
)
if USE_NUMBA:
from pyxfem.core.config import nbfloat
from pyxfem.core.nbkalgebra2d import Vec2, Sy22
import numba
import numba.core.types as nbct
import numba.typed as nbty
# pylint: disable=ungrouped-imports
from pyxfem.core.nbkalgebra2d import \
vec2_t, sy22_t, vec2_zero, sy22_zero,\
vec2_scal, vec2_square, vec2_norm, vec2_dot, vec2_add, vec2_add3, vec2_sub, \
sy22_add, sy22_mv, sy22_smv, sy22_r, sy22_r2, sy22_trace
# --- INDICES DES VARIABLES CONSERVATIVES (Tableau U) ---
DIM = 2
I_M = 0 # Masse (m)
I_MU_START = 1
I_MU_END = I_MU_START + DIM
I_ME = I_MU_END # Énergie Totale (mE)
SIZE_U = I_ME+1
# --- INDICES DES VARIABLES PRIMITIVES (Tableau W, optionnel) ---
I_RHO = 0 # Densité
I_U_START = 1 # Vitesse X
I_U_END = I_U_START + DIM # Vitesse Y
I_P = I_U_END # Pression
I_C = I_P+1 # Sound speed
SIZE_W = I_C+1
U_indexes = {"m": (I_M, I_M+1), "mu":(I_MU_START, I_MU_END), "mE":(I_ME, I_ME+1) }
U_shapes = {"size": SIZE_U, "m": (), "mu":(DIM,), "mE":() }
W_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)}
W_shapes = {"size":SIZE_W, "rho": (), "u":(DIM,), "p":(), "c":()}
@guvectorize([
(float64, float64, float64[:], float64[:,:]),
(float32, float32, float32[:], float32[:,:])
], '(),(),(m)->(m,m)', nopython=True, cache=True)
def _acoustic_impedance_tensor_2d_guvec(rho, c, n, out):
out[0,0] = rho*c*n[0]*n[0]
out[0,1] = rho*c*n[0]*n[1]
out[1,0] = rho*c*n[1]*n[0]
out[1,1] = rho*c*n[1]*n[1]
@numba.njit(parallel=True, cache = True)
def _acoustic_impedance_tensor_2d_njit(rho, c, n, out):
N = len(rho)
nv = n.shape[-2]
for j in numba.prange(N): # type ignore # pylint: disable=E1133
rhocj = rho[j] * c[j]
for r in range(nv):
njr0 = n[j,r,0]
njr1 = n[j,r,1]
out[j,r,0,0] = rhocj*njr0*njr0
out[j,r,0,1] = rhocj*njr0*njr1
out[j,r,1,0] = rhocj*njr1*njr0
out[j,r,1,1] = rhocj*njr1*njr1
[docs]
class ConservationLawManager_Euler(ConservativeLawManager):
def __init__(self, gazlaw, nbdata):
super().__init__(nbdata, U_shapes, U_indexes, W_shapes, W_indexes)
self.gazlaw = gazlaw
if USE_NUMBA:
self.numba_pressure_c_func, self.numba_soundspeed_c_func, self.fluidlaw_data = self.gazlaw.get_numba_c_funcs()
[docs]
def conservative_to_primal(self, cellmeasure):
m, mu, mE = self.get_conservative_variable(["m","mu","mE"])
rho, u, p, c = self.get_primal_variable(["rho","u", "p", "c"])
if USE_NUMBA :
pfunc = self.numba_pressure_c_func
cfunc = self.numba_soundspeed_c_func
fluidlaw_data =self.fluidlaw_data
ConservationLawManager_Euler.conservative_to_primal_numba(m, mu, mE, cellmeasure, pfunc, cfunc, fluidlaw_data, rho, u, p, c)
else:
np.divide(m, cellmeasure, out = rho)
np.divide(mu, m[..., np.newaxis], out = u)
u2 = np.square(u)
e = mE/m-fhalf*(u2[...,0] + u2[...,1])
self.gazlaw.pressure(e, rho, out = p)
self.gazlaw.soundspeed(p, rho, out = c)
if USE_NUMBA :
[docs]
@staticmethod
@numba.njit(parallel=True, cache = True)
def conservative_to_primal_numba(m, mu, mE, cell_meas, pfunc, cfunc, fluidlaw_data, rho, u, p, c):
N = len(m)
for i in numba.prange(N): # type ignore # pylint: disable=E1133
rho[i] = m[i]/cell_meas[i]
u[i, 0] = mu[i,0]/m[i]
u[i, 1] = mu[i,1]/m[i]
e = mE[i]/m[i]-fhalf * (u[i, 0]**2 + u[i, 1]**2)
p[i] = pfunc(e, rho[i], fluidlaw_data)
c[i] = cfunc(p[i], rho[i], fluidlaw_data)
[docs]
def acoustic_impedance_tensor(self, n, outZ = None):
crho, cc = self.get_primal_variable(["rho", "c"])
if outZ is None :
batch_out = np.broadcast_shapes(c.shape, n.shape[:-1])
out_dtype = n.dtype
outZ = np.empty(batch_out + (DIM, DIM), dtype=out_dtype)
# if len(n.shape[:-1]) != len(crho.shape):
# crho = crho[:, np.newaxis]
# cc = cc[:, np.newaxis]
#_acoustic_impedance_tensor_2d_guvec(crho, cc, n, outZ)
_acoustic_impedance_tensor_2d_njit(crho, cc, n, outZ)
return outZ
[docs]
def GELEES_bc_set_sliding_wall(
vids: IndexArray,
NxN: FloatArray,
ImNxN: FloatArray,
vA: FloatArray,
vb: FloatArray
) -> None:
r"""
Applies a sliding wall (slip) boundary condition to the global nodal
linear system.
This function modifies the assembled $A_r$ matrices and $b_r$ vectors
in-place for the nodes located on a sliding wall. It ensures that the
velocity vector resulting from the linear system solver will have zero
component normal to the wall ($u \cdot n = 0$).
The algebraic modification projects the system onto the tangent plane:
$A_{new} = (I - n \otimes n) A_r (I - n \otimes n) + \mathrm{Tr}(A_r)(n \otimes n)$
$b_{new} = (I - n \otimes n) b_r$
Adding the trace back to the normal component ensures the modified
matrix remains strictly invertible and well-conditioned.
Parameters
----------
vids : npt.NDArray[np.uint64]
1D array of length `(n_boundary_nodes)`. The global indices of the
vertices located on the sliding wall.
NxN : FloatArray
3D array of shape `(n_boundary_nodes, 2, 2)`. The dyadic tensor
product of the normal vectors ($n \otimes n$) at each boundary node.
ImNxN : FloatArray
3D array of shape `(n_boundary_nodes, 2, 2)`. The tangent projection
matrix ($I - n \otimes n$) at each boundary node.
vA : FloatArray
3D array of shape `(nverts, 2, 2)`. The global assembled $A_r$
matrices for all nodes in the mesh. Modified in-place.
vb : FloatArray
2D array of shape `(nverts, 2)`. The global assembled $b_r$ vectors
for all nodes in the mesh. Modified in-place.
Returns
-------
None
The global `vA` and `vb` arrays are modified in-place to reflect
the boundary constraints.
"""
boundary_vA = vA[vids]
# Extract the trace (A00 + A11) for the boundary nodes
trace_boundary_vA = boundary_vA[..., 0,0] + boundary_vA[..., 1,1]
# Apply the projection and add the penalized normal component
vA[vids] = np.matmul(np.matmul(ImNxN, boundary_vA), ImNxN) \
+ trace_boundary_vA[...,np.newaxis, np.newaxis]*NxN
# Project the RHS vector onto the tangent plane
vb[vids] = np.matvec(ImNxN, vb[vids] )
[docs]
def GELEES_bc_set_external_pressure(
vids: IndexArray,
ps: FloatArray,
mesh: GELEES_MeshStruct,
out_vA: FloatArray, # pylint: disable=w0613
out_vb: FloatArray
) -> None:
r"""Applies external pressure directly to the RHS vectors $b_{jr}$ of boundary nodes."""
for vid, p in zip(vids, ps) :
s = mesh.v2t_index[vid]
e = mesh.v2t_index[vid+1]
for cell, cv in zip(mesh.v2t_tris[s:e], mesh.v2t_v[s:e]):
out_vb[vid] -= p*mesh.C[cell, cv]
FType_v = None
FType_vc = None
if USE_NUMBA:
[docs]
@numba.njit(parallel=True, cache =True)
def rho_u_p_c_numba(U, S, pfunc, cfunc, fluidlaw_data, out_W):
N = U.shape[0]
for i in numba.prange(N): # pylint: disable=E1133
m_i, mu_i, me_i = get_conservative2d(i, U)
S_i = S[i]
rho_i = m_i/S_i
u_i = vec2_scal(1./m_i, mu_i)
e_i = me_i / m_i - npfloat(0.5) * vec2_square(u_i)
p_i = pfunc(e_i, rho_i, fluidlaw_data)
c_i = cfunc(p_i, rho_i, fluidlaw_data)
set_primal2d(rho_i, u_i, p_i, c_i, i, out_W)
# @numba.njit(cache=True, fastmath = FASTMATH)
# def GELEES_tri_dt(Sj, nCj0, nCj1, nCj2, cj):
# ''' return th critical time step in a triangle'''
# sum_normC = nCj0 + nCj1 + nCj2
# dt = Sj / sum_normC/ cj
# return dt
# @numba.njit(parallel=True, cache=True, fastmath = FASTMATH)
# def critical_timestep_numba(S, nC, c):
# min_dt = np.inf # numba does autoamatique reduction on this variable
# for j in numba.prange(len(S)): # pylint: disable=E1133
# Sj = S[j]
# nCj0 = nC[j,0]
# nCj1 = nC[j,1]
# nCj2 = nC[j,2]
# cj = c[j]
# dt = GELEES_tri_dt(Sj, nCj0, nCj1, nCj2, cj)
# # # Numba comprend automatiquement qu'il s'agit d'une réduction parallèle
# min_dt = min(min_dt, dt)
# return min_dt
[docs]
@numba.njit(parallel=False, cache=True, fastmath = FASTMATH)
def GELEES_getAbj(rhoj:npfloat, uj:Vec2, pj:npfloat, cj:npfloat, Cjr:Vec2)->Tuple[Sy22,Vec2]:
''' compute Ajr and bjr which ar part of the node k of cellj
linear systeme to solve at node jr '''
# Aj = rhoj*cj/|Cjr|*cjr X cjr
# bj = pj*Cjr + Aj.uj
scalj = rhoj * cj/vec2_norm(Cjr)
Aj = sy22_r(scalj, Cjr)
Auj = sy22_mv(Aj, uj)
bj = vec2_add( vec2_scal(pj, Cjr), Auj)
return Aj, bj
[docs]
@numba.njit(cache=True, fastmath = FASTMATH)
def bc_set_sliding_wall_numba(n, A, b):
An = sy22_mv(A, n)
scal = vec2_dot(n, An) + sy22_trace(A)
A = sy22_add(A, sy22_r(scal, n))
A = sy22_add(A, sy22_r2(npfloat(-fone), n, An))
b = vec2_add(b, vec2_scal(-vec2_dot(n, b), n))
return A, b
[docs]
def as_numba_pointer(func, signature):
"""
Encapsule une fonction Python dans une liste typée Numba.
C'est un hack nécessaire pour forcer Numba à utiliser un pointeur
générique (First-Class Function) et éviter de casser le cache JIT.
"""
# type_list_of_func = nbt.ListType(signature)
# # On instancie la liste normalement avec l'argument caché 'lsttype'
# ptr_box = nbt.List(dtype= type_list_of_func)
# ptr_box.append(func)
creator = getattr(nbty.List, "empty_list")
# 2. On instancie le conteneur physique en lui passant la signature
ptr_box = creator(signature)
# 3. Maintenant, .append() fonctionnera parfaitement !
ptr_box.append(func)
return ptr_box
data_1d_t = nbfloat[:]
# Signature pour les Nœuds : (Ar, br) = f(Ar, br, data_1d)
signature_vbc = nbct.Tuple((sy22_t, vec2_t))(sy22_t, vec2_t, data_1d_t)
# Signature pour les Coins : (Aj, bj) = f(Aj, bj, Cjr, data_1d)
signature_vcbc = nbct.Tuple((sy22_t, vec2_t))(sy22_t, vec2_t, vec2_t, data_1d_t)
# On crée les "Types" de Fonction Abstraite
FType_v = nbct.FunctionType(signature_vbc)
FType_vc = nbct.FunctionType(signature_vcbc)
[docs]
@numba.njit(signature_vbc, cache=True, fastmath=FASTMATH)
def bc_vnone_cache(Ar, br, data): # pylint: disable=W0613
return Ar, br
[docs]
@numba.njit(signature_vbc, cache=True, fastmath=FASTMATH)
def bc_sliding_wall_cache(Ar, br, data): # pylint: disable=w0613
n = (data[0], data[1])
return bc_set_sliding_wall_numba(n, Ar, br)
[docs]
@numba.njit(signature_vcbc, cache=True, fastmath=FASTMATH)
def bc_vcnone_cache(Aj, bj, Cjr, data): # pylint: disable=w0613
return Aj, bj
[docs]
@numba.njit(signature_vcbc, cache=True, fastmath=FASTMATH)
def bc_pressure_cache(Aj, bj, Cjr, data):
p = data[0]
bj_new = vec2_sub(bj, vec2_scal(p, Cjr))
return Aj, bj_new
[docs]
@numba.njit(inline='always')
def apply_bc_v(ptr_box, Ar, br, data):
""" Cache la misère du [0] pour les noeuds """
f = ptr_box[0]
return f(Ar, br, data)
[docs]
@numba.njit(inline='always')
def apply_bc_vc(ptr_box, Aj, bj, Cjr, data):
""" Cache la misère du [0] pour les coins """
f = ptr_box[0]
return f(Aj, bj, Cjr, data)
[docs]
def GELEES_bulk_vertices_loop_numba(vids, mesh, W, dt, xy, vu):
return _GELEES_bulk_vertices_loop_numba(vids, mesh, W, dt, xy, vu)
@numba.njit(parallel=True, cache=True, fastmath = FASTMATH)
def _GELEES_bulk_vertices_loop_numba(vids, mesh, W, dt, out_xy, out_vu):
# for each vertex in vids, construct the 2x2 linear system
# to find the vertex velocity and update the node position
for i in numba.prange(len(vids)): # pylint: disable=E1133
vid = vids[i]
# setup riman solver at node i
s = mesh.v2t_index[vid]
e = mesh.v2t_index[vid+1]
Ar = sy22_zero
br = vec2_zero
# Cell k connectd to node vid
for k in range(s, e):
tj = mesh.v2t_tris[k]
vr = mesh.v2t_v[k]
Cjr = (mesh.C[tj, vr, 0], mesh.C[tj, vr, 1])
rhoj, uj, pj, cj = get_primal2d(tj, W)
Aj, bj = GELEES_getAbj(rhoj, uj, pj, cj, Cjr)
Ar = sy22_add(Ar, Aj)
br = vec2_add(br, bj)
# solve rieman problem (find correct velocity)
vur = sy22_smv(Ar, br)
set_vec2(vid, out_vu, vur)
# update node position
xyr = vec2_add(get_vec2(vid, out_xy), vec2_scal(dt, vur))
set_vec2(vid, out_xy, xyr)
[docs]
def GELEES_specialized_vertices_loop(vids, mesh, W, dt, box_f_v, box_f_vc, data_v, data_vc,\
out_xy, out_vu):
_GELEES_specialized_vertices_loop(vids, mesh, W, dt, box_f_v, box_f_vc, data_v, data_vc,\
out_xy, out_vu)
@numba.njit(parallel=True, cache=True, fastmath = FASTMATH)
def _GELEES_specialized_vertices_loop(vids, mesh, W, dt, box_f_v, box_f_vc, data_v, data_vc,\
out_xy, out_vu):
# for each vertex in vids, construct the 2x2 linear system
# to find the vertex velocity and update the node position
for i in numba.prange(len(vids)): # pylint: disable=E1133
vid = vids[i]
s = mesh.v2t_index[vid]
e = mesh.v2t_index[vid+1]
Ar = sy22_zero # matrice 22 symetrique stockée ainsi (A00, A01, A11)
br = vec2_zero
node_data_v = data_v[i]
node_data_vc = data_vc[i]
# Cell k connectd to node vid
for k in range(s, e):
tj = mesh.v2t_tris[k]
vr = mesh.v2t_v[k]
Cjr = (mesh.C[tj, vr, 0], mesh.C[tj, vr, 1])
rhoj, uj, pj, cj = get_primal2d(tj, W)
Aj, bj = GELEES_getAbj(rhoj, uj, pj, cj, Cjr)
Aj, bj = apply_bc_vc(box_f_vc, Aj, bj, Cjr, node_data_vc)
Ar = sy22_add(Ar, Aj)
br = vec2_add(br, bj)
Ar, br = apply_bc_v(box_f_v, Ar, br, node_data_v)
vur = sy22_smv(Ar, br)
set_vec2(vid, out_vu, vur)
xyr = vec2_add(get_vec2(vid, out_xy), vec2_scal(dt, vur))
set_vec2(vid, out_xy, xyr)
[docs]
@numba.njit(parallel=False, cache=True, fastmath = FASTMATH)
def cell_vertex_pressure(cp, crho, cc, cu, vu, vN):
return cp - crho*cc*vec2_dot(vec2_sub(vu,cu), vN)
[docs]
def GELEES_cell_loop(dt, msh, vu, U, W, pfunc, cfunc, fluidlaw_data):
return _GELEES_cell_loop(dt, msh, vu, U, W, pfunc, cfunc, fluidlaw_data)
@numba.njit(parallel=True, cache=True, fastmath=FASTMATH)
def _GELEES_cell_loop(dt: npfloat, mesh, vu: FloatArray, U: FloatArray, W: FloatArray,\
pfunc, cfunc, fluidlaw_data) -> float:
r"""
Ultra-optimized cell loop: Fully unrolled (no inner loops),
manual prefetching to registers, minimized divisions.
"""
min_dtc = np.inf
# Extraction des pointeurs bruts au niveau global
t2v, nC, C, S, xy = mesh.t2v, mesh.nC, mesh.C, mesh.cell_meas, mesh.xy
for j in numba.prange(mesh.ncells): # pylint: disable=E1133
# =========================================================
# 1. PREFETCHING (Lecture RAM -> Registres L1)
# =========================================================
n0, n1, n2 = t2v[j, 0], t2v[j, 1], t2v[j, 2]
vur0, vur1, vur2 = get_vec2(n0, vu), get_vec2(n1, vu), get_vec2(n2, vu)
Cj0, Cj1, Cj2 = get_tri_vec2(j, C)
# Thermo courante
mj, muj, mej = get_conservative2d(j, U)
rhoj, uj, pj, cj = get_primal2d(j, W)
# =========================================================
# 2. CALCUL DES FLUX (Désenroulé, vectorisé SIMD)
# =========================================================
# Inverses des normes (3 divisions factorisées)
nC0, nC1, nC2 = get_tri_scal(j, nC)
Nj0, Nj1, Nj2 = vec2_scal(fone/nC0, Cj0),\
vec2_scal(fone/nC1, Cj1), \
vec2_scal(fone/nC2, Cj2)
# Pressions nodales (Résolveur de Riemann acoustique local)
pjr0 = cell_vertex_pressure(pj, rhoj, cj, uj, vur0, Nj0)
pjr1 = cell_vertex_pressure(pj, rhoj, cj, uj, vur1, Nj1)
pjr2 = cell_vertex_pressure(pj, rhoj, cj, uj, vur2, Nj2)
# Quantité de mouvement : q_j = sum(pjr * Cjr)
qfluxj = vec2_add3( vec2_scal(pjr0, Cj0), vec2_scal(pjr1, Cj1), vec2_scal(pjr2, Cj2))
# Énergie : e_j = sum(pjr * dot(vur, Cjr))
efluxj = pjr0*vec2_dot(vur0, Cj0) + pjr1*vec2_dot(vur1, Cj1) + pjr2*vec2_dot(vur2, Cj2)
# =========================================================
# 3. MISE A JOUR CONSERVATIVE
# =========================================================
muj_new = vec2_sub(muj, vec2_scal(dt, qfluxj) )
mej_new = mej - dt * efluxj
# =========================================================
# 4. NOUVELLE GEOMETRIE (On réutilise n0, n1, n2 !)
# =========================================================
xy0, xy1, xy2 = get_vec2(n0, xy), get_vec2(n1, xy), get_vec2(n2, xy)
Sj_new, Cj0_new, Cj1_new, Cj2_new = Cjr_triangle(xy0, xy1, xy2)
nCj0_new = vec2_norm(Cj0_new)
nCj1_new = vec2_norm(Cj1_new)
nCj2_new = vec2_norm(Cj2_new)
# =========================================================
# 5. NOUVELLE THERMODYNAMIQUE
# =========================================================
inv_Sj_new = 1.0/Sj_new
rhoj_new = mj * inv_Sj_new
inv_mj = 1.0 / mj
uj_new = vec2_scal(inv_mj, muj_new)
ej_new = mej_new * inv_mj - fhalf * vec2_square(uj_new)
pj_new = pfunc(ej_new, rhoj_new, fluidlaw_data)
cj_new = cfunc(pj_new, rhoj_new, fluidlaw_data)
# =========================================================
# 6. ECRITURE MEMOIRE ET PAS DE TEMPS
# =========================================================
set_conservative2d(mj, muj_new, mej_new, j, U)
set_primal2d(rhoj_new, uj_new, pj_new, cj_new, j, W)
S[j] = Sj_new
set_tri_scal(j, nC, nCj0_new, nCj1_new, nCj2_new)
set_tri_vec2(j, C, Cj0_new, Cj1_new, Cj2_new)
# Pas de temps critic
dtc = Sj_new / ((nCj0_new + nCj1_new + nCj2_new) * cj_new)
min_dtc = min(min_dtc, dtc)
return min_dtc
[docs]
@numba.njit(inline='always', fastmath = FASTMATH)
def get_vec2(index, array):
return (array[index, 0], array[index, 1])
[docs]
@numba.njit(cache=True, fastmath = FASTMATH)
def get_tri_vec2(index, array):
return (array[index, 0, 0], array[index,0, 1]), \
(array[index, 1, 0], array[index,1, 1]), \
(array[index, 2, 0], array[index,2, 1]),
[docs]
@numba.njit(inline='always', fastmath = FASTMATH)
def get_tri_scal(index, array):
return array[index, 0], array[index, 1], array[index, 2]
[docs]
@numba.njit(inline='always', fastmath = FASTMATH)
def set_vec2(index, array, v):
array[index, 0], array[index, 1] = v[0], v[1]
[docs]
@numba.njit(inline='always', fastmath = FASTMATH)
def set_tri_scal(index, array, s0, s1, s2):
array[index, 0], array[index, 1], array[index, 2] = s0, s1, s2
[docs]
@numba.njit(cache= True, fastmath = FASTMATH)
def set_tri_vec2(index, array, v0, v1, v2):
array[index, 0, 0], array[index, 0, 1] = v0[0], v0[1]
array[index, 1, 0], array[index, 1, 1] = v1[0], v1[1]
array[index, 2, 0], array[index, 2, 1] = v2[0], v2[1]
[docs]
@numba.njit(cache=True, fastmath = FASTMATH)
def get_trixy(elemid, mesh):
tri = mesh.t2v[elemid]
n0, n1, n2 = tri[0], tri[1], tri[2]
xy = mesh.xy
xy0 = get_vec2(n0, xy)
xy1 = get_vec2(n1, xy)
xy2 = get_vec2(n2, xy)
return xy0, xy1, xy2
[docs]
@numba.njit(cache=True, fastmath = FASTMATH)
def get_conservative2d(index, U):
m, mu, me = U[index, I_M], (U[index, I_MU_START], U[index, I_MU_START+1]), U[index, I_ME]
return m, mu, me
[docs]
@numba.njit(cache=True, fastmath = FASTMATH)
def set_conservative2d(m, mu, me, index, U):
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
[docs]
@numba.njit(cache=True, fastmath = FASTMATH)
def get_primal2d( index, W):
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]
return rho, u, p, c
[docs]
@numba.njit(cache=True, fastmath = FASTMATH)
def set_primal2d(rho, u, p, c, index, W):
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
[docs]
class GlaceExplicitLagrangian_EulerEquationSolver(GlaceExplicitLagrangianBaseSolver):# pylint: disable=C0103
"""
A 2D Explicit Lagrangian Finite Volume Solver for the Euler Equations
using the GLACE (Godunov-type Lagrangian Scheme for Compressible Euler) scheme.
This class manages the mesh topology, physical state variables, boundary
conditions, and time integration (both in pure NumPy and Numba JIT modes).
"""
def __init__(
self,
xy: FloatArray,
tris: IndexArray,
boundary_conditions: Dict[str, Any],
use_numba: bool = True
) -> None:
"""
Initializes the GLACE solver.
Parameters
----------
xy : FloatArray
2D array of shape `(nverts, 2)` containing the initial nodal coordinates.
tris : IndexArray
2D array of shape `(ncells, 3)` containing the cell-to-vertex connectivity.
boundary_conditions : Dict[str, Any]
Dictionary containing the boundary condition definitions (types, nodes,
normals, and prescribed pressures).
use_numba : bool, optional
If True, the solver will attempt to use the optimized Numba JIT-compiled
routines. Defaults to True.
"""
self.use_numba = False
if use_numba :
if USE_NUMBA:
self.use_numba = True
else :
print("Warning : Trying to use numba will the __FILE__ USE_NUMBA is set to false")
self.gaz = PerfectGaz()
mesh = setup_MeshStruct(xy, tris)
nverts = mesh.nverts
ncells = mesh.ncells
conservation_law = ConservationLawManager_Euler(self.gaz, ncells)
super().__init__(conservation_law, mesh)
self.t = 0.
self.vu = np.zeros((nverts, DIM), dtype=npfloat)
# conservative vars (m, mux, muy, me) per cell
# primary vars (rho, ux, uy, p, c) per cell
self._init_bc(boundary_conditions)
if self.use_numba :
#computeCjr_triangles_numba(self.mesh)
self.numba_pressure_c_func, self.numba_soundspeed_c_func, self.fluidlaw_data = self.gaz.get_numba_c_funcs()
self.step = self.step_numba
else:
#Cjr_triangle_np(self.mesh)
self.step = self.step_np
def _init_bc(self, boundary_conditions: Dict[str, Any]) -> None:
# 1. ACCUMULATION : On récolte les infos physiques par noeud
# Chaque noeud a un dictionnaire avec ses normales et ses pressions
node_data = defaultdict(lambda: {'normals': [], 'pressures': []})
for bc_data in boundary_conditions.values():
bc_type = bc_data["type"]
vids = bc_data["vid"]
if bc_type == "sliding_wall":
n = np.array(bc_data["n"], dtype=npfloat)
for vid in vids:
node_data[vid]['normals'].append(n)
elif bc_type == "external_pressure":
p = float(bc_data["pressure"])
for vid in vids:
node_data[vid]['pressures'].append(p)
# 2. CLASSIFICATION : Préparation des listes par catégorie
categories = {
"sliding_wall": {"nodes": [], "N": []},
"external_pressure": {"nodes": [], "pressures": []},
"pressure_and_wall": {"nodes": [], "pressures": [], "N": []},
"fixed_point": {"nodes": []}
}
for vid, data in node_data.items():
normals = data['normals']
pressures = data['pressures']
n_walls = len(normals)
n_press = len(pressures)
if n_walls == 1 and n_press == 0:
categories["sliding_wall"]["nodes"].append(vid)
categories["sliding_wall"]["N"].append(normals[0])
elif n_walls == 0 and n_press >= 1:
categories["external_pressure"]["nodes"].append(vid)
categories["external_pressure"]["pressures"].append(np.mean(pressures))
elif n_walls == 1 and n_press >= 1:
categories["pressure_and_wall"]["nodes"].append(vid)
categories["pressure_and_wall"]["N"].append(normals[0])
categories["pressure_and_wall"]["pressures"].append(np.mean(pressures))
elif n_walls >= 2:
# 2 murs ou plus qui se croisent = point fixe (coin)
categories["fixed_point"]["nodes"].append(vid)
else:
print(f"Warning: BC non gérée pour le noeud {vid}")
# 3. VECTORISATION : Création du dictionnaire final avec les arrays NumPy
self.bc = {}
for cat_name, cat_data in categories.items():
bc_dict ={} #{"nodes": np.array(cat_data["nodes"], dtype=IndexType)}
bc_dict["nodes"] = np.array(cat_data["nodes"], dtype=IndexType)
if "N" in cat_data:
N_arr = np.array(cat_data["N"], dtype=npfloat).reshape(-1, 2)
bc_dict["N"] = N_arr
# Pour le solveur NumPy pur, on vectorise la création des tenseurs (NxN et I-NxN)
if not self.use_numba:
if len(N_arr) > 0:
NxN = N_arr[:, :, np.newaxis] * N_arr[:, np.newaxis, :]
bc_dict["NxN"] = NxN
bc_dict["ImNxN"] = np.eye(2) - NxN
else:
# Ensures the keys exist even if there are no nodes
bc_dict["NxN"] = np.empty((0, 2, 2), dtype=npfloat)
bc_dict["ImNxN"] = np.empty((0, 2, 2), dtype=npfloat)
if "pressures" in cat_data:
bc_dict["pressures"] = np.array(cat_data["pressures"], dtype=npfloat)
self.bc[cat_name] = bc_dict
# Separating Bulk from boundary vertices
boundary_vids = np.concatenate([bc["nodes"] for bc in self.bc.values()])
self.bulk_vids = np.setdiff1d(np.arange(self.mesh.nverts), boundary_vids)
# BINDING NUMBA : attaching function pointer and data for each bc type.
if self.use_numba:
for bc_type, bc_dict in self.bc.items():
n_nodes = len(bc_dict["nodes"])
if bc_type == "sliding_wall":
bc_dict["vdata_cache"] = bc_dict["N"]
bc_dict["vfun_cache"] = as_numba_pointer(bc_sliding_wall_cache, FType_v)
bc_dict["vcdata_cache"] = np.empty((n_nodes, 0), dtype=npfloat)
bc_dict["vcfun_cache"] = as_numba_pointer(bc_vcnone_cache, FType_vc)
elif bc_type == "external_pressure":
bc_dict["vdata_cache"] = np.empty((n_nodes, 0), dtype=npfloat)
bc_dict["vfun_cache"] = as_numba_pointer(bc_vnone_cache, FType_v)
bc_dict["vcdata_cache"] = bc_dict.get("pressures", np.empty(0, dtype=npfloat)).reshape(-1, 1)
bc_dict["vcfun_cache"] = as_numba_pointer(bc_pressure_cache, FType_vc)
elif bc_type == "pressure_and_wall":
bc_dict["vdata_cache"] = bc_dict["N"]
bc_dict["vfun_cache"] = as_numba_pointer(bc_sliding_wall_cache, FType_v)
bc_dict["vcdata_cache"] = bc_dict.get("pressures", np.empty(0, dtype=npfloat)).reshape(-1, 1)
bc_dict["vcfun_cache"] = as_numba_pointer(bc_pressure_cache, FType_vc)
# Nothing done for "fixed_point" : we just have to make sure those points are never moved.
def _set_boundary_conditions(self)-> None:
# external_pressure
bc_data = self.bc["external_pressure"]
vids = bc_data["nodes"]
ps = bc_data["pressures"]
GELEES_bc_set_external_pressure(vids, ps, self.mesh, self.vA, self.vb)
# sliding wall :
bc_data = self.bc["sliding_wall"]
vids = bc_data["nodes"]
NxN = bc_data["NxN"]
ImNxN = bc_data["ImNxN"]
GELEES_bc_set_sliding_wall(vids, NxN, ImNxN, self.vA, self.vb)
# pressure_and_wall
bc_data = self.bc["pressure_and_wall"]
vids = bc_data["nodes"]
ps = bc_data["pressures"]
NxN = bc_data["NxN"]
ImNxN = bc_data["ImNxN"]
GELEES_bc_set_external_pressure(vids, ps, self.mesh, self.vA, self.vb)
GELEES_bc_set_sliding_wall(vids, NxN, ImNxN, self.vA, self.vb)
# Fixed (2 sliding walls)
bc_data = self.bc["fixed_point"]
vids = bc_data["nodes"]
self.vb[vids] = 0.
def _riemann_solver_setup(self)-> None:
cvn: FloatArray = self.mesh.N
cvnC: FloatArray = self.mesh.nC
cvC: FloatArray = self.mesh.C
cvA: FloatArray = self.cvA
cvb: FloatArray = self.cvb
cu, cp = self.conservation_law.get_primal_variable(["u", "p"])
self.conservation_law.acoustic_impedance_tensor(cvn, outZ = self.cvA)
if USE_NUMBA :
GlaceExplicitLagrangian_EulerEquationSolver._riemann_solver_setup_njit(cu, cp, cvA, cvC, cvnC, cvA, cvb)
else:
cvA *= cvnC[..., np.newaxis, np.newaxis]
cp_ = np.expand_dims(cp, axis = (1,2))
# b_jr = p * Cjr + A_jr * u
cvb[:] = cvC*cp_ \
+ self.cvA[...,0]*cu[:, np.newaxis, np.newaxis, 0] \
+ self.cvA[...,1]*cu[:, np.newaxis, np.newaxis, 1]
if USE_NUMBA :
@staticmethod
@numba.njit(parallel=True, cache = True, fastmath = FASTMATH)
def _riemann_solver_setup_njit(cu, cp, cvZ, cvC, cvnC, out_cvA, out_cvb):
N =cu.shape[0]
for j in numba.prange(N): # type ignore # pylint: disable=E1133
pj = cp[j]
uj0 = cu[j,0]
uj1 = cu[j,1]
vAj = out_cvA[j]
vZj = cvZ[j]
vnCj = cvnC[j]
vbj = out_cvb[j]
vCj = cvC[j]
for r in range(3):
Cjr = vCj[r]
nCjr = vnCj[r]
Zjr = vZj[r]
Ajr = vAj[r]
Ajr[0, 0] = nCjr*Zjr[0 ,0]
Ajr[0, 1] = nCjr*Zjr[0 ,1]
Ajr[1, 0] = nCjr*Zjr[1 ,0]
Ajr[1, 1] = nCjr*Zjr[1 ,1]
bjr = vbj[r]
bjr[0] = pj*Cjr[0] + Ajr[0, 0]*uj0 + Ajr[0, 1]*uj1
bjr[1] = pj*Cjr[1] + Ajr[1, 0]*uj0 + Ajr[1, 1]*uj1
[docs]
def set_initial_state(
self,
density: Callable[[FloatArray], FloatArray],
velocity: Callable[[FloatArray], FloatArray],
pressure: Callable[[FloatArray], FloatArray]
) -> None:
"""
Sets the initial macroscopic state of the fluid for all cells.
Parameters
----------
density : Callable[[FloatArray], FloatArray]
A function taking cell center coordinates (N, 2) and returning
an array of initial cu, cp = self.conservation_law.get_primal_variable(["u", "p"])
cp_ = np.expand_dims(cp, axis = (1,2))
# b_jr = p * Cjr + A_jr * u
self.cvb[:] = cvC*cp_ \
+ self.cvA[...,0]*cu[:, np.newaxis, np.newaxis, 0] \
+ self.cvA[...,1]*cu[:, np.newaxis, np.newaxis, 1]
densities (N,).
velocity : Callable[[FloatArray], FloatArray]
A function taking cell center coordinates (N, 2) and returning
an array of initial velocity vectors (N, 2).
pressure : Callable[[FloatArray], FloatArray]
A function taking cell center coordinates (N, 2) and returning
an array of initial pressures (N,).
"""
xytris = self.mesh.xy[self.mesh.t2v]
xycog = np.sum(xytris, axis =1)/3.
densitycog = np.array(density(xycog), dtype = npfloat)
velocitycog = np.array(velocity(xycog), dtype = npfloat)
pressurecog = np.array(pressure(xycog), dtype = npfloat)
cell_measure = self.mesh.cell_meas
m = densitycog*cell_measure
int_spec_enrg = self.gaz.internal_specific_energy(pressurecog,densitycog)
kin_spec_enrg = fhalf*(velocitycog[:, 0]**2 + velocitycog[:, 1]**2)
getU = self.conservation_law.get_conservative_variable
getU("m")[:] = m
getU("mu")[:] = velocitycog*m[..., np.newaxis]
getU("mE")[:] = (int_spec_enrg+kin_spec_enrg)*m
self.conservation_law.conservative_to_primal(cell_measure)
[docs]
def critical_timestep(self) -> float:
"""
Computes the maximum allowable critical time step for the entire mesh
based on the CFL condition (acoustic wave speeds).
Returns
-------
dtmax : float
The critical time step (dt).
"""
V = self.mesh.cell_meas
cc = self.get_primal_variable("c") # per cell sound speed
nC = self.mesh.nC
return type(self).critical_timestep_pnjit(V, nC, cc)
# if True : #self.use_numba:
# nC = self.mesh.nC
# return critical_timestep_numba(V, nC, cc)
# else:
# C = self.mesh.C
# cvnormC = np.sqrt( np.square(C[...,0]) + np.square(C[...,1]) )
# #sum_cvnormC =np.sum(cvnormC, axis =1)
# sum_cvnormC = cvnormC[:,0] + cvnormC[:,1] + cvnormC[:,2]
# cL = V/sum_cvnormC
# dtmax = np.min(cL/cc)
# return dtmax
[docs]
def step(self, dt: float) -> float: # pylint: disable=E0202:method-hidden
"""
Advances the simulation by one time step `dt`.
Notes
-----
This method is dynamically routed at runtime. Depending on the `use_numba`
flag provided during the solver's initialization, calling `step()` will
execute either:
- `step_numba()` : The heavily optimized JIT-compiled solver.
- `step_np()` : The pure vectorized NumPy solver.
Parameters
----------
dt : float
The duration of the time step.
Returns
-------
dtnew : float
The new computed critical time step for the next iteration.
"""
# This line will never be executed thanks to dynamic affctation to self.step in __init__
raise NotImplementedError("This method is dynamically assigned during __init__.")
[docs]
def step_np(self, dt: float) -> float:
"""
Advances the simulation by one time step `dt` using the pure NumPy solver.
Parameters
----------
dt : float
The duration of the time step.
Returns
-------
dtcnew : float
The new computed critical time step for the next iteration.
"""
#nc = self.ncells
vu, cvqflux = self._riemann_solver()
mesh = self.mesh
cvu = self.cvu
#graph.extract_txy(vu, self.mesh.t2v, cvu)
#cvu = vu[self.mesh.t2v]
GlaceExplicitLagrangian_EulerEquationSolver.update_conservative(cvqflux, cvu, dt, self.conservation_law)
#moving node
mesh.xy[:] += dt*vu
# update primal variable (on moved cells)
update_MeshStruct(mesh.xy, mesh)
self.conservation_law.conservative_to_primal(mesh.cell_meas)
#update time
self.t += dt
#compute critical time step
dtcnew = self.critical_timestep()
return dtcnew
[docs]
@staticmethod
def update_conservative(cvqflux, cvu, dt, cons_law_man):
mu, mE = cons_law_man.get_conservative_variable(["mu", "mE"])
if USE_NUMBA:
GlaceExplicitLagrangian_EulerEquationSolver.update_conservative_njit(cvqflux, cvu, dt, mu, mE)
else:
mu[:, 0] -= dt*(cvqflux[:,0,0] + cvqflux[:,1,0] + cvqflux[:,2,0])
mu[:, 1] -= dt*(cvqflux[:,0,1] + cvqflux[:,1,1] + cvqflux[:,2,1])
mE -= dt*( cvqflux[:,0,0] * cvu[:,0,0] + cvqflux[:,0,1] * cvu[:,0,1] \
+ cvqflux[:,1,0] * cvu[:,1,0] + cvqflux[:,1,1] * cvu[:,1,1] \
+ cvqflux[:,2,0] * cvu[:,2,0] + cvqflux[:,2,1] * cvu[:,2,1])
if USE_NUMBA:
[docs]
@staticmethod
@numba.njit(parallel = True, cache = True, fastmath = FASTMATH )
def update_conservative_njit(cvqflux, cvu, dt, mu, mE):
for j in numba.prange(cvqflux.shape[0]): # type ignore # pylint: disable=E1133
vqfluxj = cvqflux[j]
vuj = cvu[j]
mu[j, 0] -= dt*(vqfluxj[0,0] + vqfluxj[1,0] + vqfluxj[2,0])
mu[j, 1] -= dt*(vqfluxj[0,1] + vqfluxj[1,1] + vqfluxj[2,1])
mE[j] -= dt*( vqfluxj[0,0] * vuj[0,0] + vqfluxj[0,1] * vuj[0,1] \
+ vqfluxj[1,0] * vuj[1,0] + vqfluxj[1,1] * vuj[1,1] \
+ vqfluxj[2,0] * vuj[2,0] + vqfluxj[2,1] * vuj[2,1])
if USE_NUMBA :
[docs]
def step_numba(self, dt: float) -> float:
"""
Advances the simulation by one time step `dt` using the optimized
Numba JIT-compiled solver.
Parameters
----------
dt : float
The duration of the time step.
Returns
-------
dtnew : float
The new computed critical time step for the next iteration.
"""
mesh = self.mesh # the msh
W = self.conservation_law.W # cell primal variable
U = self.conservation_law.U # cell conservative variable
vu = self.vu # vertex velocity
# 1. bulk nodes update
bulk_vids = self.bulk_vids
#GELEES_vertices_loop(bulk_vids, mesh, W, dt, bc_vnone, bc_vcnone, mesh.xy, vu)
GELEES_bulk_vertices_loop_numba(bulk_vids, mesh, W, dt, mesh.xy, vu)
# 2. boundary node update
# 2.1 wall bc
for bc_type in ["sliding_wall", "external_pressure", "pressure_and_wall"]:
bc = self.bc[bc_type]
bc_vids = bc["nodes"]
if len(bc_vids) > 0 :
box_bc_vfun = bc["vfun_cache"]
box_bc_vcfun = bc["vcfun_cache"]
v_data = bc["vdata_cache"]
vc_data = bc["vcdata_cache"]
GELEES_specialized_vertices_loop(bc_vids, mesh, W, dt, box_bc_vfun, box_bc_vcfun, v_data, vc_data, mesh.xy, vu)
p_c_func = self.numba_pressure_c_func
sp_c_func = self.numba_soundspeed_c_func
fluidlaw_data = self.fluidlaw_data
dtnew = GELEES_cell_loop(dt, mesh, vu, U, W, p_c_func, sp_c_func, fluidlaw_data)
self.t += dt
return dtnew