Coverage for pyxfem/solidmechanics/finite_volume_2d_hyperelasticity.py: 0%

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

8 

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

13 

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 typing import Any, Callable, Dict, Tuple, NamedTuple, Sequence, List, Union, assert_never, Optional, overload 

21from collections import defaultdict 

22from numba.core.bytecode import get_code_object 

23import numpy as np 

24 

25from numba import guvectorize, float32, float64 

26import numba 

27 

28from pyxfem.core.graph import T2V, extract_txy 

29from pyxfem.core.config import USE_NUMBA, FASTMATH, IndexType, npfloat, fzero, fone, fhalf,\ 

30 FloatArray, IndexArray 

31from pyxfem.core.smallsize_tensor import gemm22_kernel_njit, max_eigval22sym, gemm22, inv, sysmv22, AAT22_kernel, max_eigval22sym_kernel_njit, axpy 

32from pyxfem.core.finite_volume_2d import t_tangentmap, GELEES_MeshStruct, setup_MeshStruct, update_MeshStruct, ConservativeLawManager, \ 

33 GlaceExplicitLagrangianBaseSolver 

34 

35 

36from pyxfem.solidmechanics.largestrain import acoustic_impedance_tensor, leftCauchyGreen, acoustic_impedance_tensor_pnjit #wavespeed, acoustic_tensor_space 

37from pyxfem.solidmechanics.hyperelasticity import HyperElasticMaterial 

38 

39 

40DIM = 2 

41# --- INDICES DES VARIABLES CONSERVATIVES (Tableau U) --- 

42I_M = 0 # Masse (m) 

43I_MU_START = 1 # impulsion : mass x velocity vector 

44I_MU_END = I_MU_START+DIM 

45I_ME = I_MU_END # Énergie Totale (mE) 

46I_MF_START = I_ME +1 

47I_MF_END = I_MF_START+DIM*DIM # cell mass multiplied by current position gradient F = dx/dX 

48SIZE_U = I_MF_END 

49 

50# --- INDICES DES VARIABLES PRIMITIVES (Tableau W) --- 

51I_RHO = 0 # Densité 

52I_U_START = 1 

53I_U_END = I_U_START + DIM # Cauchy Stress xx 

54I_SIG_START = I_U_END #Cauchy Stress xx 

55I_SIG_END = I_SIG_START + DIM*DIM#_END Cauchy Stress 

56I_C_START = I_SIG_END# starting index for the spacial tangent operator 

57I_C_END = I_C_START + DIM*DIM*DIM*DIM 

58SIZE_W = I_C_END 

59 

60U_indexes = {"m": (I_M, I_M+1), "mu":(I_MU_START, I_MU_END), "mE":(I_ME, I_ME+1), "mF":(I_MF_START, I_MF_END) } 

61U_shapes = {"size": SIZE_U, "m": (), "mu":(DIM,), "mE":(), "mF":(DIM, DIM) } 

62 

63W_indexes = {"rho": (I_RHO, I_RHO+1), "u": (I_U_START, I_U_END), "sigma": (I_SIG_START, I_SIG_END), "c":(I_C_START, I_C_END)} 

64W_shapes = {"size":SIZE_W, "rho": (), "u":(DIM,), "sigma":(DIM,DIM), "c":(DIM, DIM, DIM, DIM)} 

65 

66class ConservationLawManager_Hyperelasticity(ConservativeLawManager): 

67 def __init__(self, hyplaw, nbdata): 

68 super().__init__(nbdata, U_shapes, U_indexes, W_shapes, W_indexes) 

69 self.hyplaw = hyplaw 

70 

71 def conservative_to_primal(self, cell_measure): 

72 """ Compute from U (the conservative variable) the primal Variable W depending on the law""" 

73 m, mu, mF = self.get_conservative_variable(["m", "mu","mF"]) 

74 rho, u, sigma, c = self.get_primal_variable(["rho", "u", "sigma", "c"]) 

75 if USE_NUMBA: 

76 props_sigma_c = self.hyplaw.get_properties() 

77 fun_sigma_c = type(self.hyplaw).sigma_c_kernel 

78 ConservationLawManager_Hyperelasticity.conservative_to_primal_(m, mu, mF, cell_measure, fun_sigma_c, props_sigma_c, rho, u, sigma, c) 

79 else : 

80 F = sigma # using F spce to compute sigma 

81 np.divide(m, cell_measure, out = rho) 

82 np.divide(mu, m[..., np.newaxis], out = u) 

83 np.divide(mF, m[..., np.newaxis,np.newaxis], F) 

84 b = leftCauchyGreen(F) 

85 self.hyplaw.Cauchy(b, out_sigma = sigma) 

86 self.hyplaw.spacial_elasticty_tensor(b, out_c = c) 

87 

88 @staticmethod 

89 @numba.njit(parallel = True, cache =True, fastmath = FASTMATH) 

90 def conservative_to_primal_(m, mu, mF, cell_measure, fun_sigma_c, props_sigma_c, rho, u, sigma, c): 

91 n = m.shape[0] 

92 F = sigma 

93 for j in numba.prange(n): 

94 rho[j] = m[j]/cell_measure[j] 

95 u[j] = mu[j]/m[j] 

96 F[j] = mF[j]/m[j] 

97 fun_sigma_c(F[j], props_sigma_c, sigma[j], c[j]) 

98 

99 def acoustic_impedance_tensor(self, n, outZ = None): 

100 """ Compute the impedance """ 

101 crho, csig, cc = self.get_primal_variable(["rho","sigma","c"]) 

102 cc = cc[:,np.newaxis] 

103 csig = csig[:,np.newaxis] 

104 crho = crho[:, np.newaxis] 

105 if outZ is None: 

106 batch_out = np.broadcast_shapes(cc.shape[:-4], csig.shape[:-2], n.shape[:-1]) 

107 out_dtype = n.dtype 

108 outZ = np.empty(batch_out + (2, 2), dtype=out_dtype) 

109 return acoustic_impedance_tensor_pnjit(crho, cc, csig, n, outZ) 

110 #return acoustic_impedance_tensor(crho, cc, csig, n, outZ) 

111 

112 

113@numba.njit(cache = True) 

114def update_conservative_vars_kernel(cvqflux, cvu, dt, in_out_U_cell): 

115 q_x = cvqflux[0, 0] + cvqflux[1, 0] + cvqflux[2, 0] 

116 q_y = cvqflux[0, 1] + cvqflux[1, 1] + cvqflux[2, 1] 

117 e_flux = (cvqflux[0, 0] * cvu[0, 0] + cvqflux[0, 1] * cvu[0, 1]) + \ 

118 (cvqflux[1, 0] * cvu[1, 0] + cvqflux[1, 1] * cvu[1, 1]) + \ 

119 (cvqflux[2, 0] * cvu[2, 0] + cvqflux[2, 1] * cvu[2, 1]) 

120 in_out_U_cell[I_MU_START] -= dt * q_x 

121 in_out_U_cell[I_MU_START+1] -= dt * q_y 

122 in_out_U_cell[I_ME] -= dt * e_flux 

123 

124@numba.njit(parallel=True, cache = True) 

125def update_conservative_vars_pnjit(cvqflux, cvu, dt, in_out_U): 

126 for j in numba.prange(cvqflux.shape[0]): 

127 update_conservative_vars_kernel(cvqflux[j], cvu[j], dt, in_out_U[j]) 

128 

129@guvectorize([ 

130 (float64[:,:], float64[:,:], float64, float64[:]), # 5 types 

131 (float32[:,:], float32[:,:], float32, float32[:]) 

132], '(n_v,dim),(n_v,dim),(), (n_u)', nopython=True, cache=True) 

133def update_conservative_vars_guvec(cvqflux, cvu, dt, in_out_U_cell): 

134 ''' cvqflux : cell vertex flux (stress) 

135 cvu : cell vertex velocity 

136 U_cell : cell conservative Variables 

137 Note this verdion do not update mF the defomation multiplyed by cell mass 

138 ''' 

139 # flux computation 

140 update_conservative_vars_kernel(cvqflux, cvu, dt, in_out_U_cell) 

141 

142 

143 

144 

145 

146 

147class Additional_MeshStruct(NamedTuple): 

148 XY: FloatArray 

149 tXY: FloatArray 

150 F0: FloatArray 

151 txy: FloatArray 

152 iF0: FloatArray 

153 F1: FloatArray 

154 F: FloatArray 

155 

156 

157def setup_AdditionalMeshStruct(XY, xy, tris): 

158 tXY = extract_txy(XY, tris) 

159 F0 = t_tangentmap(tXY) 

160 iF0 = inv(F0) 

161 txy = extract_txy(xy, tris) 

162 F1 = t_tangentmap(txy) 

163 F = gemm22(F1, iF0) 

164 return Additional_MeshStruct( XY= XY, tXY = tXY, F0 = F0, iF0 = iF0, txy = txy, F1 =F1, F = F) 

165 

166def update_AdditionalMeshStruct(mesh: GELEES_MeshStruct, additional_mesh_data: Additional_MeshStruct): 

167 amesh = additional_mesh_data 

168 if USE_NUMBA: 

169 update_AdditionalMeshStruct_pnjit(mesh, amesh) 

170 else: 

171 extract_txy(mesh.xy, mesh.t2v, out = amesh.txy) 

172 t_tangentmap(amesh.txy, out = amesh.F1) 

173 gemm22( amesh.F1, amesh.iF0, amesh.F) 

174 

175@numba.njit(parallel = True, cache =True, fastmath = True) 

176def update_AdditionalMeshStruct_pnjit(mesh: GELEES_MeshStruct, additional_mesh_data: Additional_MeshStruct): 

177 amesh = additional_mesh_data 

178 xy = mesh.xy 

179 nc = mesh.ncells 

180 for j in numba.prange(nc): # type ignore # pylint: disable=E1133 

181 v0, v1, v2 = mesh.t2v[j] 

182 xy0 = xy[v0] 

183 xy1 = xy[v1] 

184 xy2 = xy[v2] 

185 invF0 = amesh.iF0[j] 

186 F = amesh.F[j] 

187 #F1 = amesh.F1[j] 

188 F1 = F # using F as a temporary storage for F1 

189 F1[0,0] = xy1[0] - xy0[0] 

190 F1[1,0] = xy1[1] - xy0[1] 

191 F1[0,1] = xy2[0] - xy0[0] 

192 F1[1,1] = xy2[1] - xy0[1] 

193 gemm22_kernel_njit(F1, invF0, F) 

194 

195 

196 

197 

198 

199 

200 

201 

202 

203 

204 

205def GELEES_bc_set_sliding_wall( 

206 vids: IndexArray, 

207 NxN: FloatArray, 

208 ImNxN: FloatArray, 

209 vA: FloatArray, 

210 vb: FloatArray 

211) -> None: 

212 r""" 

213 Applies a sliding wall (slip) boundary condition to the global nodal 

214 linear system. 

215 

216 This function modifies the assembled $A_r$ matrices and $b_r$ vectors 

217 in-place for the nodes located on a sliding wall. It ensures that the 

218 velocity vector resulting from the linear system solver will have zero 

219 component normal to the wall ($u \cdot n = 0$). 

220 

221 The algebraic modification projects the system onto the tangent plane: 

222 $A_{new} = (I - n \otimes n) A_r (I - n \otimes n) + \mathrm{Tr}(A_r)(n \otimes n)$ 

223 $b_{new} = (I - n \otimes n) b_r$ 

224 

225 Adding the trace back to the normal component ensures the modified 

226 matrix remains strictly invertible and well-conditioned. 

227 

228 Parameters 

229 ---------- 

230 vids : npt.NDArray[np.uint64] 

231 1D array of length `(n_boundary_nodes)`. The global indices of the 

232 vertices located on the sliding wall. 

233 NxN : FloatArray 

234 3D array of shape `(n_boundary_nodes, 2, 2)`. The dyadic tensor 

235 product of the normal vectors ($n \otimes n$) at each boundary node. 

236 ImNxN : FloatArray 

237 3D array of shape `(n_boundary_nodes, 2, 2)`. The tangent projection 

238 matrix ($I - n \otimes n$) at each boundary node. 

239 vA : FloatArray 

240 3D array of shape `(nverts, 2, 2)`. The global assembled $A_r$ 

241 matrices for all nodes in the mesh. Modified in-place. 

242 vb : FloatArray 

243 2D array of shape `(nverts, 2)`. The global assembled $b_r$ vectors 

244 for all nodes in the mesh. Modified in-place. 

245 

246 Returns 

247 ------- 

248 None 

249 The global `vA` and `vb` arrays are modified in-place to reflect 

250 the boundary constraints. 

251 """ 

252 boundary_vA = vA[vids] 

253 # Extract the trace (A00 + A11) for the boundary nodes 

254 trace_boundary_vA = boundary_vA[..., 0,0] + boundary_vA[..., 1,1] 

255 # Apply the projection and add the penalized normal component 

256 vA[vids] = np.matmul(np.matmul(ImNxN, boundary_vA), ImNxN) \ 

257 + trace_boundary_vA[...,np.newaxis, np.newaxis]*NxN 

258 # Project the RHS vector onto the tangent plane 

259 vb[vids] = np.matvec(ImNxN, vb[vids] ) 

260 

261def GELEES_bc_set_external_pressure( 

262 vids: IndexArray, 

263 ps: FloatArray, 

264 mesh: GELEES_MeshStruct, 

265 out_vA: FloatArray, # pylint: disable=w0613 

266 out_vb: FloatArray 

267) -> None: 

268 r"""Applies external pressure directly to the RHS vectors $b_{jr}$ of boundary nodes.""" 

269 for vid, p in zip(vids, ps) : 

270 s = mesh.v2t_index[vid] 

271 e = mesh.v2t_index[vid+1] 

272 for cell, cv in zip(mesh.v2t_tris[s:e], mesh.v2t_v[s:e]): 

273 out_vb[vid] -= p*mesh.C[cell, cv] 

274 

275FType_v = None 

276FType_vc = None 

277 

278# Signature de dimension : 

279# (), (), () -> scalaires (lamb, mu, rho) 

280# (r, d), (r) -> géométrie : normales aux noeuds C_jr et leurs normes (3, 2) et (3,) 

281# (d, d), (d) -> variables de l'écoulement : contrainte sigma (2, 2) et vitesse u (2,) 

282# (r, d, d), (r, d)-> sorties modifiées sur place : cvA (3, 2, 2) et cvb (3, 2) 

283@guvectorize( 

284 [(float64, float64, float64, float64[:, :], float64[:], float64[:, :], float64[:], float64[:, :, :], float64[:, :])], 

285 '(),(),(),(r,d),(r),(d,d),(d)->(r,d,d),(r,d)', 

286 nopython=True, 

287 target='parallel' 

288) 

289def riemann_setup_1_guvec(lamb, mu, crho, cvC, cvnC, csig, cu, cvA, cvb): 

290 rho_aL = np.sqrt((lamb + 2.0 * mu) * crho) 

291 rho_aT = np.sqrt(mu * crho) 

292 

293 # Nombre de noeuds par cellule (généralement 3) 

294 r_len = cvC.shape[0] 

295 

296 for r in range(r_len): 

297 C0 = cvC[r, 0] 

298 C1 = cvC[r, 1] 

299 nC = cvnC[r] 

300 

301 # Factorisation de la division pour éviter les coûts CPU 

302 # On ajoute une sécurité numérique contre la division par zéro 

303 inv_nC = 1.0 / nC if nC > 1e-15 else 0.0 

304 

305 C00 = C0 * C0 * inv_nC 

306 C01 = C0 * C1 * inv_nC 

307 C11 = C1 * C1 * inv_nC 

308 

309 # --- 1. Calcul de la matrice A_jr --- 

310 cvA[r, 0, 0] = rho_aL * C00 + rho_aT * (nC - C00) 

311 cvA[r, 0, 1] = rho_aL * C01 - rho_aT * C01 

312 cvA[r, 1, 0] = cvA[r, 0, 1] # Exploitation de la symétrie 

313 cvA[r, 1, 1] = rho_aL * C11 + rho_aT * (nC - C11) 

314 

315 # --- 2. Calcul du vecteur b_jr --- 

316 # Équivalent de : -(csig @ Cjr) 

317 sigC_0 = csig[0, 0] * C0 + csig[0, 1] * C1 

318 sigC_1 = csig[1, 0] * C0 + csig[1, 1] * C1 

319 

320 # Équivalent de : A_jr @ u_j 

321 Au_0 = cvA[r, 0, 0] * cu[0] + cvA[r, 0, 1] * cu[1] 

322 Au_1 = cvA[r, 1, 0] * cu[0] + cvA[r, 1, 1] * cu[1] 

323 

324 # Résultat final b_jr = -sig . C + A . u 

325 cvb[r, 0] = -sigC_0 + Au_0 

326 cvb[r, 1] = -sigC_1 + Au_1 

327 

328class GlaceExplicitLagrangian_HyperElasticty_EquationSolver(GlaceExplicitLagrangianBaseSolver): # pylint: disable=C0103 

329 """ 

330 A 2D Explicit Lagrangian Finite Volume Solver for the Euler Equations 

331 using the GLACE (Godunov-type Lagrangian Scheme for Compressible Euler) scheme. 

332 

333 This class manages the mesh topology, physical state variables, boundary 

334 conditions, and time integration (both in pure NumPy and Numba JIT modes). 

335 """ 

336 def __init__( 

337 self, 

338 law: HyperElasticMaterial, 

339 XY: FloatArray, 

340 xy: FloatArray, 

341 tris: IndexArray, 

342 boundary_conditions: Dict[str, Any], 

343 use_numba: bool = False 

344 ) -> None: 

345 """ 

346 Initializes the GLACE solver. 

347 

348 Parameters 

349 ---------- 

350 xy : FloatArray 

351 2D array of shape `(nverts, 2)` containing the initial nodal coordinates. 

352 tris : IndexArray 

353 2D array of shape `(ncells, 3)` containing the cell-to-vertex connectivity. 

354 boundary_conditions : Dict[str, Any] 

355 Dictionary containing the boundary condition definitions (types, nodes, 

356 normals, and prescribed pressures). 

357 use_numba : bool, optional 

358 If True, the solver will attempt to use the optimized Numba JIT-compiled 

359 routines. Defaults to True. 

360 """ 

361 

362 self.law = law 

363 ne = len(tris) 

364 conservation_law = ConservationLawManager_Hyperelasticity(law,ne) 

365 super().__init__(conservation_law, mesh = setup_MeshStruct(xy, tris)) 

366 

367 self.mesh_adata = setup_AdditionalMeshStruct(XY, xy, tris) 

368 ncells = self.mesh.ncells 

369 nverts = self.mesh.nverts 

370 

371 self.t = 0. 

372 self.vu = np.zeros((nverts, DIM), dtype=npfloat) 

373 # boundary condition setup 

374 self._init_bc(boundary_conditions) 

375 self.cvZ = np.zeros((ncells, 3, 2, 2), dtype=npfloat) #Impedance tensors 

376 self.cc = np.zeros((ncells), dtype=npfloat) 

377 self.step = self.step_np 

378 

379 def _init_bc(self, boundary_conditions: Dict[str, Any]) -> None: 

380 # 1. ACCUMULATION : On récolte les infos physiques par noeud 

381 # Chaque noeud a un dictionnaire avec ses normales et ses pressions 

382 node_data = defaultdict(lambda: {'normals': [], 'pressures': []}) 

383 for bc_data in boundary_conditions.values(): 

384 bc_type = bc_data["type"] 

385 vids = bc_data["vid"] 

386 

387 if bc_type == "sliding_wall": 

388 n = np.array(bc_data["n"], dtype=npfloat) 

389 for vid in vids: 

390 node_data[vid]['normals'].append(n) 

391 elif bc_type == "external_pressure": 

392 p = float(bc_data["pressure"]) 

393 for vid in vids: 

394 node_data[vid]['pressures'].append(p) 

395 # 2. CLASSIFICATION : Préparation des listes par catégorie 

396 categories = { 

397 "sliding_wall": {"nodes": [], "N": []}, 

398 "external_pressure": {"nodes": [], "pressures": []}, 

399 "pressure_and_wall": {"nodes": [], "pressures": [], "N": []}, 

400 "fixed_point": {"nodes": []} 

401 } 

402 

403 for vid, data in node_data.items(): 

404 normals = data['normals'] 

405 pressures = data['pressures'] 

406 n_walls = len(normals) 

407 n_press = len(pressures) 

408 if n_walls == 1 and n_press == 0: 

409 categories["sliding_wall"]["nodes"].append(vid) 

410 categories["sliding_wall"]["N"].append(normals[0]) 

411 elif n_walls == 0 and n_press >= 1: 

412 categories["external_pressure"]["nodes"].append(vid) 

413 categories["external_pressure"]["pressures"].append(np.mean(pressures)) 

414 elif n_walls == 1 and n_press >= 1: 

415 categories["pressure_and_wall"]["nodes"].append(vid) 

416 categories["pressure_and_wall"]["N"].append(normals[0]) 

417 categories["pressure_and_wall"]["pressures"].append(np.mean(pressures)) 

418 elif n_walls >= 2: 

419 # 2 murs ou plus qui se croisent = point fixe (coin) 

420 categories["fixed_point"]["nodes"].append(vid) 

421 else: 

422 print(f"Warning: BC non gérée pour le noeud {vid}") 

423 

424 # 3. VECTORISATION : Création du dictionnaire final avec les arrays NumPy 

425 self.bc = {} 

426 for cat_name, cat_data in categories.items(): 

427 bc_dict ={} #{"nodes": np.array(cat_data["nodes"], dtype=IndexType)} 

428 bc_dict["nodes"] = np.array(cat_data["nodes"], dtype=IndexType) 

429 if "N" in cat_data: 

430 N_arr = np.array(cat_data["N"], dtype=npfloat).reshape(-1, 2) 

431 bc_dict["N"] = N_arr 

432 # Pour le solveur NumPy pur, on vectorise la création des tenseurs (NxN et I-NxN) 

433 # if not self.use_numba: 

434 # if len(N_arr) > 0: 

435 # NxN = N_arr[:, :, np.newaxis] * N_arr[:, np.newaxis, :] 

436 # bc_dict["NxN"] = NxN 

437 # bc_dict["ImNxN"] = np.eye(2) - NxN 

438 # else: 

439 # Ensures the keys exist even if there are no nodes 

440 bc_dict["NxN"] = np.empty((0, 2, 2), dtype=npfloat) 

441 bc_dict["ImNxN"] = np.empty((0, 2, 2), dtype=npfloat) 

442 

443 if "pressures" in cat_data: 

444 bc_dict["pressures"] = np.array(cat_data["pressures"], dtype=npfloat) 

445 self.bc[cat_name] = bc_dict 

446 # Separating Bulk from boundary vertices 

447 boundary_vids = np.concatenate([bc["nodes"] for bc in self.bc.values()]) 

448 self.bulk_vids = np.setdiff1d(np.arange(self.mesh.nverts), boundary_vids) 

449 # BINDING NUMBA : attaching function pointer and data for each bc type. 

450 # if self.use_numba: 

451 # for bc_type, bc_dict in self.bc.items(): 

452 # n_nodes = len(bc_dict["nodes"]) 

453 # if bc_type == "sliding_wall": 

454 # bc_dict["vdata_cache"] = bc_dict["N"] 

455 # bc_dict["vfun_cache"] = as_numba_pointer(bc_sliding_wall_cache, FType_v) 

456 # bc_dict["vcdata_cache"] = np.empty((n_nodes, 0), dtype=npfloat) 

457 # bc_dict["vcfun_cache"] = as_numba_pointer(bc_vcnone_cache, FType_vc) 

458 # elif bc_type == "external_pressure": 

459 # bc_dict["vdata_cache"] = np.empty((n_nodes, 0), dtype=npfloat) 

460 # bc_dict["vfun_cache"] = as_numba_pointer(bc_vnone_cache, FType_v) 

461 # bc_dict["vcdata_cache"] = bc_dict.get("pressures", np.empty(0, dtype=npfloat)).reshape(-1, 1) 

462 # bc_dict["vcfun_cache"] = as_numba_pointer(bc_pressure_cache, FType_vc) 

463 # elif bc_type == "pressure_and_wall": 

464 # bc_dict["vdata_cache"] = bc_dict["N"] 

465 # bc_dict["vfun_cache"] = as_numba_pointer(bc_sliding_wall_cache, FType_v) 

466 # bc_dict["vcdata_cache"] = bc_dict.get("pressures", np.empty(0, dtype=npfloat)).reshape(-1, 1) 

467 # bc_dict["vcfun_cache"] = as_numba_pointer(bc_pressure_cache, FType_vc) 

468 # # Nothing done for "fixed_point" : we just have to make sure those points are never moved. 

469 

470 def _set_boundary_conditions(self)-> None: 

471 # external_pressure 

472 bc_data = self.bc["external_pressure"] 

473 vids = bc_data["nodes"] 

474 ps = bc_data["pressures"] 

475 GELEES_bc_set_external_pressure(vids, ps, self.mesh, self.vA, self.vb) 

476 # sliding wall : 

477 bc_data = self.bc["sliding_wall"] 

478 vids = bc_data["nodes"] 

479 NxN = bc_data["NxN"] 

480 ImNxN = bc_data["ImNxN"] 

481 GELEES_bc_set_sliding_wall(vids, NxN, ImNxN, self.vA, self.vb) 

482 # pressure_and_wall 

483 bc_data = self.bc["pressure_and_wall"] 

484 vids = bc_data["nodes"] 

485 ps = bc_data["pressures"] 

486 NxN = bc_data["NxN"] 

487 ImNxN = bc_data["ImNxN"] 

488 GELEES_bc_set_external_pressure(vids, ps, self.mesh, self.vA, self.vb) 

489 GELEES_bc_set_sliding_wall(vids, NxN, ImNxN, self.vA, self.vb) 

490 # Fixed (2 sliding walls) 

491 bc_data = self.bc["fixed_point"] 

492 vids = bc_data["nodes"] 

493 self.vb[vids] = 0. 

494 

495 def compute_impedance_tensors(self, mesh:GELEES_MeshStruct, out_cvZ:FloatArray|None = None): 

496 """ Compute the impdence tensor Z, per cell vertex (cv)""" 

497 cvn : FloatArray = self.mesh.N 

498 return self.conservation_law.acoustic_impedance_tensor(cvn, out_cvZ) 

499 

500 @staticmethod 

501 def cell_soundspeed(cvZ, cc): 

502 if USE_NUMBA: 

503 GlaceExplicitLagrangian_HyperElasticty_EquationSolver.cell_soundspeed_pnjit(cvZ, cc) 

504 else: 

505 cc[:] = np.amax(max_eigval22sym(cvZ), axis = -1) 

506 

507 @staticmethod 

508 @numba.njit(parallel=True, cache=True,fastmath = FASTMATH) 

509 def cell_soundspeed_pnjit(cvZ, cc): 

510 for j in numba.prange(cvZ.shape[0]): 

511 vZj = cvZ[j] 

512 cj = 0. 

513 for r in range(cvZ.shape[1]): 

514 cj = max(max_eigval22sym_kernel_njit(vZj[r]), cj) 

515 cc[j] = cj 

516 

517 def set_initial_state( 

518 self, 

519 #density: Callable[[FloatArray], FloatArray], 

520 velocity: Callable[[FloatArray], FloatArray] 

521 #pressure: Callable[[FloatArray], FloatArray] 

522 ) -> None: 

523 """ 

524 Sets the initial macroscopic state of the fluid for all cells. 

525 

526 Parameters 

527 ---------- 

528 density : Callable[[FloatArray], FloatArray] 

529 A function taking cell center coordinates (N, 2) and returning 

530 an array of initial densities (N,). 

531 velocity : Callable[[FloatArray], FloatArray] 

532 A function taking cell center coordinates (N, 2) and returning 

533 an array of initial velocity vectors (N, 2). 

534 pressure : Callable[[FloatArray], FloatArray] 

535 A function taking cell center coordinates (N, 2) and returning 

536 an array of initial pressures (N,). 

537 """ 

538 xytris = self.mesh.xy[self.mesh.t2v] 

539 xycog = np.sum(xytris, axis =1)/3. 

540 #densitycog = np.array(density(xycog), dtype = npfloat) 

541 velocitycog = np.array(velocity(xycog), dtype = npfloat) 

542 #pressurecog = np.array(pressure(xycog), dtype = npfloat) 

543 

544 cell_measure = self.mesh.cell_meas 

545 m = self.law.rho*cell_measure 

546 getU = self.conservation_law.get_conservative_variable 

547 getU("m")[:] = m 

548 getU("mu")[:] = velocitycog*m[..., np.newaxis] 

549 kin_spec_enrg = fhalf*(velocitycog[:, 0]**2 + velocitycog[:, 1]**2) 

550 getU("mE")[:] = kin_spec_enrg*m 

551 getU("mF")[:] = self.mesh_adata.F*m[:,np.newaxis, np.newaxis] 

552 self.conservation_law.conservative_to_primal(cell_measure) 

553 self.cvZ = self.compute_impedance_tensors(self.mesh) 

554 self.cc[:] = np.amax(max_eigval22sym(self.cvZ), axis = -1) 

555 

556 def _riemann_solver_setup_0(self)-> None: 

557 r""" 

558 Using simpliyed GLACE ... 

559 Computes the local cell-vertex matrix and vector contributions for the 

560 nodal acoustic Riemann solver. 

561 

562 This function vectorizes the calculation of the $A_{jr}$ matrices and 

563 $b_{jr}$ vectors for all cells and their vertices simultaneously. These 

564 local contributions will later be assembled into a global linear system 

565 at each node to compute the nodal velocities. 

566 

567 The continuous formulations for the acoustic solver are: 

568 $A_{jr} = \rho_j c_j \frac{C_{jr} \otimes C_{jr}}{|C_{jr}|}$ 

569 $b_{jr} = p_j C_{jr} + A_{jr} \cdot u_j$ 

570 

571 Parameters 

572 ---------- 

573 cvC : FloatArray 

574 3D array of shape `(ncells, 3, 2)` containing the corner normal vectors 

575 $C_{jr}$ for each vertex of each cell. 

576 cvnC : FloatArray 

577 2D array of shape `(ncells, 3)` containing the Euclidean norms (magnitudes) 

578 of the corner normal vectors $|C_{jr}|$. 

579 W : FloatArray 

580 2D array of shape `(ncells, 5)` containing the primitive variables for 

581 each cell (density, x-velocity, y-velocity, pressure, sound speed). 

582 outcvA : FloatArray 

583 4D output array of shape `(ncells, 3, 2, 2)`. The computed symmetric 

584 $2 \times 2$ matrices $A_{jr}$ are stored here. Modified in-place. 

585 outcvb : FloatArray 

586 3D output array of shape `(ncells, 3, 2)`. The computed 2D vectors 

587 $b_{jr}$ are stored here. Modified in-place. 

588 

589 Returns 

590 ------- 

591 None 

592 The outputs `outcvA` and `outcvb` are modified in-place to avoid memory 

593 reallocation.""" 

594 # GELEES_riemann_solver_setup( self.mesh.C, self.mesh.nC, self.W, self.cvA, self.cvb) 

595 cvC: FloatArray = self.mesh.C 

596 cvnC: FloatArray = self.mesh.nC 

597 #cvZ: FloatArray = self.cvZ 

598 csig: FloatArray = self.get_primal_variable("sigma") 

599 cu : FloatArray = self.get_primal_variable("u") 

600 

601 crho = self.get_primal_variable("rho") 

602 lamb = self.law.lamb 

603 mu = self.law.mu 

604 

605 for j in range(self.cvA.shape[0]): 

606 rhoj = crho[j] 

607 #print(rhoj) 

608 aLj = np.sqrt((lamb+2*mu)/rhoj) 

609 aTj = np.sqrt(mu/rhoj) 

610 vCj = cvC[j] 

611 vnCj = cvnC[j] 

612 for r in range(self.cvA.shape[1]): 

613 Cjr = vCj[r] 

614 nCjr = vnCj[r] 

615 CjrxCjr_00 = Cjr[0]*Cjr[0] 

616 CjrxCjr_01 = Cjr[0]*Cjr[1] 

617 CjrxCjr_10 = CjrxCjr_01 

618 CjrxCjr_11 = Cjr[1]*Cjr[1] 

619 self.cvA[j,r][0,0] = rhoj *aLj * CjrxCjr_00/nCjr +rhoj*aTj*( nCjr - CjrxCjr_00/nCjr ) 

620 self.cvA[j,r][0,1] = rhoj *aLj * CjrxCjr_01/nCjr - rhoj*aTj*( CjrxCjr_01/nCjr ) 

621 self.cvA[j,r][1,0] = rhoj *aLj * CjrxCjr_10/nCjr - rhoj*aTj*( CjrxCjr_10/nCjr ) 

622 self.cvA[j,r][1,1] = rhoj *aLj * CjrxCjr_11/nCjr +rhoj*aTj*( nCjr - CjrxCjr_11/nCjr ) 

623 

624 #self.cvA[:]= cvZ*cvnC[..., np.newaxis, np.newaxis] 

625 

626 # b_jr = sigj . Cjr + A_jr . u 

627 csig = csig[:,np.newaxis,...] 

628 cvC = cvC[..., np.newaxis] 

629 self.cvb[:] = -(csig@cvC).squeeze() 

630 

631 self.cvb[:] += self.cvA[...,0]*cu[:, np.newaxis, np.newaxis, 0] \ 

632 + self.cvA[...,1]*cu[:, np.newaxis, np.newaxis, 1] 

633 

634 def _riemann_solver_setup_1(self) -> None : 

635 

636 lamb = self.law.lamb 

637 mu = self.law.mu 

638 cvC: FloatArray = self.mesh.C 

639 cvnC: FloatArray = self.mesh.nC 

640 crho, csig, cu = self.get_primal_variable(["rho", "sigma", "u"]) 

641 cvA, cvb = self.cvA, self.cvb 

642 riemann_setup_1_guvec(lamb, mu, crho, cvC, cvnC, csig, cu, cvA, cvb) 

643 

644 

645 def _riemann_solver_setup(self)-> None: 

646 r""" 

647 Computes the local cell-vertex matrix and vector contributions for the 

648 nodal acoustic Rie 

649mann solver. 

650 

651 This function vectorizes the calculation of the $A_{jr}$ matrices and 

652 $b_{jr}$ vectors for all cells and their vertices simultaneously. These 

653 local contributions will later be assembled into a global linear system 

654 at each node to compute the nodal velocities. 

655 

656 The continuous formulations for the acoustic solver are: 

657 $A_{jr} = \rho_j c_j \frac{C_{jr} \otimes C_{jr}}{|C_{jr}|}$ 

658 $b_{jr} = p_j C_{jr} + A_{jr} \cdot u_j$ 

659 

660 Parameters 

661 ---------- 

662 cvC : FloatArray 

663 3D array of shape `(ncells, 3, 2)` containing the corner normal vectors 

664 $C_{jr}$ for each vertex of each cell. 

665 cvnC : FloatArray 

666 2D array of shape `(ncells, 3)` containing the Euclidean norms (magnitudes) 

667 of the corner normal vectors $|C_{jr}|$. 

668 W : FloatArray 

669 2D array of shape `(ncells, 5)` containing the primitive variables for 

670 each cell (density, x-velocity, y-velocity, pressure, sound speed). 

671 outcvA : FloatArray 

672 4D output array of shape `(ncells, 3, 2, 2)`. The computed symmetric 

673 $2 \times 2$ matrices $A_{jr}$ are stored here. Modified in-place. 

674 outcvb : FloatArray 

675 3D output array of shape `(ncells, 3, 2)`. The computed 2D vectors 

676 $b_{jr}$ are stored here. Modified in-place. 

677 

678 Returns 

679 ------- 

680 None 

681 The outputs `outcvA` and `outcvb` are modified in-place to avoid memory 

682 reallocation.""" 

683 # GELEES_riemann_solver_setup( self.mesh.C, self.mesh.nC, self.W, self.cvA, self.cvb) 

684 cvC: FloatArray = self.mesh.C 

685 cvnC: FloatArray = self.mesh.nC 

686 cvZ: FloatArray = self.cvZ 

687 cvA: FloatArray = self.cvA 

688 cvb: FloatArray = self.cvb 

689 csig: FloatArray = self.get_primal_variable("sigma") 

690 cu : FloatArray = self.get_primal_variable("u") 

691 if USE_NUMBA: 

692 GlaceExplicitLagrangian_HyperElasticty_EquationSolver._riemann_solver_setup_njit(cu, csig, cvZ, cvC, cvnC, cvA, cvb) 

693 else: 

694 cvA[:]= cvZ*cvnC[..., np.newaxis, np.newaxis] 

695 # b_jr = sigj . Cjr + A_jr . u 

696 csig = csig[:,np.newaxis,...] 

697 cvC = cvC[..., np.newaxis] 

698 cvb[:] = -(csig@cvC).squeeze() 

699 cvb[:] += self.cvA[...,0]*cu[:, np.newaxis, np.newaxis, 0] \ 

700 + self.cvA[...,1]*cu[:, np.newaxis, np.newaxis, 1] 

701 if USE_NUMBA : 

702 @staticmethod 

703 @numba.njit(parallel=True, cache = True, fastmath = FASTMATH) 

704 def _riemann_solver_setup_njit(cu, csig, cvZ, cvC, cvnC, out_cvA, out_cvb): 

705 N =cu.shape[0] 

706 for j in numba.prange(N): 

707 sigj = csig[j] 

708 uj0 = cu[j,0] 

709 uj1 = cu[j,1] 

710 vAj = out_cvA[j] 

711 vZj = cvZ[j] 

712 vnCj = cvnC[j]