#!/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 typing import Any, Callable, Dict, Tuple, NamedTuple, Sequence, List, Union, assert_never, Optional, overload
from collections import defaultdict
from numba.core.bytecode import get_code_object
import numpy as np
from numba import guvectorize, float32, float64
import numba
from pyxfem.core.graph import T2V, extract_txy
from pyxfem.core.config import USE_NUMBA, FASTMATH, IndexType, npfloat, fzero, fone, fhalf,\
FloatArray, IndexArray
from pyxfem.core.smallsize_tensor import gemm22_kernel_njit, max_eigval22sym, gemm22, inv, sysmv22, AAT22_kernel, max_eigval22sym_kernel_njit, axpy
from pyxfem.core.finite_volume_2d import t_tangentmap, GELEES_MeshStruct, setup_MeshStruct, update_MeshStruct, ConservativeLawManager, \
GlaceExplicitLagrangianBaseSolver
from pyxfem.solidmechanics.largestrain import acoustic_impedance_tensor, leftCauchyGreen, acoustic_impedance_tensor_pnjit #wavespeed, acoustic_tensor_space
from pyxfem.solidmechanics.hyperelasticity import HyperElasticMaterial
DIM = 2
# --- INDICES DES VARIABLES CONSERVATIVES (Tableau U) ---
I_M = 0 # Masse (m)
I_MU_START = 1 # impulsion : mass x velocity vector
I_MU_END = I_MU_START+DIM
I_ME = I_MU_END # Énergie Totale (mE)
I_MF_START = I_ME +1
I_MF_END = I_MF_START+DIM*DIM # cell mass multiplied by current position gradient F = dx/dX
SIZE_U = I_MF_END
# --- INDICES DES VARIABLES PRIMITIVES (Tableau W) ---
I_RHO = 0 # Densité
I_U_START = 1
I_U_END = I_U_START + DIM # Cauchy Stress xx
I_SIG_START = I_U_END #Cauchy Stress xx
I_SIG_END = I_SIG_START + DIM*DIM#_END Cauchy Stress
I_C_START = I_SIG_END# starting index for the spacial tangent operator
I_C_END = I_C_START + DIM*DIM*DIM*DIM
SIZE_W = I_C_END
U_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) }
U_shapes = {"size": SIZE_U, "m": (), "mu":(DIM,), "mE":(), "mF":(DIM, DIM) }
W_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)}
W_shapes = {"size":SIZE_W, "rho": (), "u":(DIM,), "sigma":(DIM,DIM), "c":(DIM, DIM, DIM, DIM)}
[docs]
class ConservationLawManager_Hyperelasticity(ConservativeLawManager):
def __init__(self, hyplaw, nbdata):
super().__init__(nbdata, U_shapes, U_indexes, W_shapes, W_indexes)
self.hyplaw = hyplaw
[docs]
def conservative_to_primal(self, cell_measure):
""" Compute from U (the conservative variable) the primal Variable W depending on the law"""
m, mu, mF = self.get_conservative_variable(["m", "mu","mF"])
rho, u, sigma, c = self.get_primal_variable(["rho", "u", "sigma", "c"])
if USE_NUMBA:
props_sigma_c = self.hyplaw.get_properties()
fun_sigma_c = type(self.hyplaw).sigma_c_kernel
ConservationLawManager_Hyperelasticity.conservative_to_primal_(m, mu, mF, cell_measure, fun_sigma_c, props_sigma_c, rho, u, sigma, c)
else :
F = sigma # using F spce to compute sigma
np.divide(m, cell_measure, out = rho)
np.divide(mu, m[..., np.newaxis], out = u)
np.divide(mF, m[..., np.newaxis,np.newaxis], F)
b = leftCauchyGreen(F)
self.hyplaw.Cauchy(b, out_sigma = sigma)
self.hyplaw.spacial_elasticty_tensor(b, out_c = c)
[docs]
@staticmethod
@numba.njit(parallel = True, cache =True, fastmath = FASTMATH)
def conservative_to_primal_(m, mu, mF, cell_measure, fun_sigma_c, props_sigma_c, rho, u, sigma, c):
n = m.shape[0]
F = sigma
for j in numba.prange(n):
rho[j] = m[j]/cell_measure[j]
u[j] = mu[j]/m[j]
F[j] = mF[j]/m[j]
fun_sigma_c(F[j], props_sigma_c, sigma[j], c[j])
[docs]
def acoustic_impedance_tensor(self, n, outZ = None):
""" Compute the impedance """
crho, csig, cc = self.get_primal_variable(["rho","sigma","c"])
cc = cc[:,np.newaxis]
csig = csig[:,np.newaxis]
crho = crho[:, np.newaxis]
if outZ is None:
batch_out = np.broadcast_shapes(cc.shape[:-4], csig.shape[:-2], n.shape[:-1])
out_dtype = n.dtype
outZ = np.empty(batch_out + (2, 2), dtype=out_dtype)
return acoustic_impedance_tensor_pnjit(crho, cc, csig, n, outZ)
#return acoustic_impedance_tensor(crho, cc, csig, n, outZ)
[docs]
@numba.njit(cache = True)
def update_conservative_vars_kernel(cvqflux, cvu, dt, in_out_U_cell):
q_x = cvqflux[0, 0] + cvqflux[1, 0] + cvqflux[2, 0]
q_y = cvqflux[0, 1] + cvqflux[1, 1] + cvqflux[2, 1]
e_flux = (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])
in_out_U_cell[I_MU_START] -= dt * q_x
in_out_U_cell[I_MU_START+1] -= dt * q_y
in_out_U_cell[I_ME] -= dt * e_flux
[docs]
@numba.njit(parallel=True, cache = True)
def update_conservative_vars_pnjit(cvqflux, cvu, dt, in_out_U):
for j in numba.prange(cvqflux.shape[0]):
update_conservative_vars_kernel(cvqflux[j], cvu[j], dt, in_out_U[j])
@guvectorize([
(float64[:,:], float64[:,:], float64, float64[:]), # 5 types
(float32[:,:], float32[:,:], float32, float32[:])
], '(n_v,dim),(n_v,dim),(), (n_u)', nopython=True, cache=True)
def update_conservative_vars_guvec(cvqflux, cvu, dt, in_out_U_cell):
''' cvqflux : cell vertex flux (stress)
cvu : cell vertex velocity
U_cell : cell conservative Variables
Note this verdion do not update mF the defomation multiplyed by cell mass
'''
# flux computation
update_conservative_vars_kernel(cvqflux, cvu, dt, in_out_U_cell)
[docs]
class Additional_MeshStruct(NamedTuple):
XY: FloatArray
tXY: FloatArray
F0: FloatArray
txy: FloatArray
iF0: FloatArray
F1: FloatArray
F: FloatArray
[docs]
def setup_AdditionalMeshStruct(XY, xy, tris):
tXY = extract_txy(XY, tris)
F0 = t_tangentmap(tXY)
iF0 = inv(F0)
txy = extract_txy(xy, tris)
F1 = t_tangentmap(txy)
F = gemm22(F1, iF0)
return Additional_MeshStruct( XY= XY, tXY = tXY, F0 = F0, iF0 = iF0, txy = txy, F1 =F1, F = F)
[docs]
def update_AdditionalMeshStruct(mesh: GELEES_MeshStruct, additional_mesh_data: Additional_MeshStruct):
amesh = additional_mesh_data
if USE_NUMBA:
update_AdditionalMeshStruct_pnjit(mesh, amesh)
else:
extract_txy(mesh.xy, mesh.t2v, out = amesh.txy)
t_tangentmap(amesh.txy, out = amesh.F1)
gemm22( amesh.F1, amesh.iF0, amesh.F)
[docs]
@numba.njit(parallel = True, cache =True, fastmath = True)
def update_AdditionalMeshStruct_pnjit(mesh: GELEES_MeshStruct, additional_mesh_data: Additional_MeshStruct):
amesh = additional_mesh_data
xy = mesh.xy
nc = mesh.ncells
for j in numba.prange(nc): # type ignore # pylint: disable=E1133
v0, v1, v2 = mesh.t2v[j]
xy0 = xy[v0]
xy1 = xy[v1]
xy2 = xy[v2]
invF0 = amesh.iF0[j]
F = amesh.F[j]
#F1 = amesh.F1[j]
F1 = F # using F as a temporary storage for F1
F1[0,0] = xy1[0] - xy0[0]
F1[1,0] = xy1[1] - xy0[1]
F1[0,1] = xy2[0] - xy0[0]
F1[1,1] = xy2[1] - xy0[1]
gemm22_kernel_njit(F1, invF0, F)
[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
# Signature de dimension :
# (), (), () -> scalaires (lamb, mu, rho)
# (r, d), (r) -> géométrie : normales aux noeuds C_jr et leurs normes (3, 2) et (3,)
# (d, d), (d) -> variables de l'écoulement : contrainte sigma (2, 2) et vitesse u (2,)
# (r, d, d), (r, d)-> sorties modifiées sur place : cvA (3, 2, 2) et cvb (3, 2)
@guvectorize(
[(float64, float64, float64, float64[:, :], float64[:], float64[:, :], float64[:], float64[:, :, :], float64[:, :])],
'(),(),(),(r,d),(r),(d,d),(d)->(r,d,d),(r,d)',
nopython=True,
target='parallel'
)
def riemann_setup_1_guvec(lamb, mu, crho, cvC, cvnC, csig, cu, cvA, cvb):
rho_aL = np.sqrt((lamb + 2.0 * mu) * crho)
rho_aT = np.sqrt(mu * crho)
# Nombre de noeuds par cellule (généralement 3)
r_len = cvC.shape[0]
for r in range(r_len):
C0 = cvC[r, 0]
C1 = cvC[r, 1]
nC = cvnC[r]
# Factorisation de la division pour éviter les coûts CPU
# On ajoute une sécurité numérique contre la division par zéro
inv_nC = 1.0 / nC if nC > 1e-15 else 0.0
C00 = C0 * C0 * inv_nC
C01 = C0 * C1 * inv_nC
C11 = C1 * C1 * inv_nC
# --- 1. Calcul de la matrice A_jr ---
cvA[r, 0, 0] = rho_aL * C00 + rho_aT * (nC - C00)
cvA[r, 0, 1] = rho_aL * C01 - rho_aT * C01
cvA[r, 1, 0] = cvA[r, 0, 1] # Exploitation de la symétrie
cvA[r, 1, 1] = rho_aL * C11 + rho_aT * (nC - C11)
# --- 2. Calcul du vecteur b_jr ---
# Équivalent de : -(csig @ Cjr)
sigC_0 = csig[0, 0] * C0 + csig[0, 1] * C1
sigC_1 = csig[1, 0] * C0 + csig[1, 1] * C1
# Équivalent de : A_jr @ u_j
Au_0 = cvA[r, 0, 0] * cu[0] + cvA[r, 0, 1] * cu[1]
Au_1 = cvA[r, 1, 0] * cu[0] + cvA[r, 1, 1] * cu[1]
# Résultat final b_jr = -sig . C + A . u
cvb[r, 0] = -sigC_0 + Au_0
cvb[r, 1] = -sigC_1 + Au_1
[docs]
class GlaceExplicitLagrangian_HyperElasticty_EquationSolver(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,
law: HyperElasticMaterial,
XY: FloatArray,
xy: FloatArray,
tris: IndexArray,
boundary_conditions: Dict[str, Any],
use_numba: bool = False
) -> 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.law = law
ne = len(tris)
conservation_law = ConservationLawManager_Hyperelasticity(law,ne)
super().__init__(conservation_law, mesh = setup_MeshStruct(xy, tris))
self.mesh_adata = setup_AdditionalMeshStruct(XY, xy, tris)
ncells = self.mesh.ncells
nverts = self.mesh.nverts
self.t = 0.
self.vu = np.zeros((nverts, DIM), dtype=npfloat)
# boundary condition setup
self._init_bc(boundary_conditions)
self.cvZ = np.zeros((ncells, 3, 2, 2), dtype=npfloat) #Impedance tensors
self.cc = np.zeros((ncells), dtype=npfloat)
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.
[docs]
def compute_impedance_tensors(self, mesh:GELEES_MeshStruct, out_cvZ:FloatArray|None = None):
""" Compute the impdence tensor Z, per cell vertex (cv)"""
cvn : FloatArray = self.mesh.N
return self.conservation_law.acoustic_impedance_tensor(cvn, out_cvZ)
[docs]
@staticmethod
def cell_soundspeed(cvZ, cc):
if USE_NUMBA:
GlaceExplicitLagrangian_HyperElasticty_EquationSolver.cell_soundspeed_pnjit(cvZ, cc)
else:
cc[:] = np.amax(max_eigval22sym(cvZ), axis = -1)
[docs]
@staticmethod
@numba.njit(parallel=True, cache=True,fastmath = FASTMATH)
def cell_soundspeed_pnjit(cvZ, cc):
for j in numba.prange(cvZ.shape[0]):
vZj = cvZ[j]
cj = 0.
for r in range(cvZ.shape[1]):
cj = max(max_eigval22sym_kernel_njit(vZj[r]), cj)
cc[j] = cj
[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 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 = self.law.rho*cell_measure
getU = self.conservation_law.get_conservative_variable
getU("m")[:] = m
getU("mu")[:] = velocitycog*m[..., np.newaxis]
kin_spec_enrg = fhalf*(velocitycog[:, 0]**2 + velocitycog[:, 1]**2)
getU("mE")[:] = kin_spec_enrg*m
getU("mF")[:] = self.mesh_adata.F*m[:,np.newaxis, np.newaxis]
self.conservation_law.conservative_to_primal(cell_measure)
self.cvZ = self.compute_impedance_tensors(self.mesh)
self.cc[:] = np.amax(max_eigval22sym(self.cvZ), axis = -1)
def _riemann_solver_setup_0(self)-> None:
r"""
Using simpliyed GLACE ...
Computes the local cell-vertex matrix and vector contributions for the
nodal acoustic Riemann solver.
This function vectorizes the calculation of the $A_{jr}$ matrices and
$b_{jr}$ vectors for all cells and their vertices simultaneously. These
local contributions will later be assembled into a global linear system
at each node to compute the nodal velocities.
The continuous formulations for the acoustic solver are:
$A_{jr} = \rho_j c_j \frac{C_{jr} \otimes C_{jr}}{|C_{jr}|}$
$b_{jr} = p_j C_{jr} + A_{jr} \cdot u_j$
Parameters
----------
cvC : FloatArray
3D array of shape `(ncells, 3, 2)` containing the corner normal vectors
$C_{jr}$ for each vertex of each cell.
cvnC : FloatArray
2D array of shape `(ncells, 3)` containing the Euclidean norms (magnitudes)
of the corner normal vectors $|C_{jr}|$.
W : FloatArray
2D array of shape `(ncells, 5)` containing the primitive variables for
each cell (density, x-velocity, y-velocity, pressure, sound speed).
outcvA : FloatArray
4D output array of shape `(ncells, 3, 2, 2)`. The computed symmetric
$2 \times 2$ matrices $A_{jr}$ are stored here. Modified in-place.
outcvb : FloatArray
3D output array of shape `(ncells, 3, 2)`. The computed 2D vectors
$b_{jr}$ are stored here. Modified in-place.
Returns
-------
None
The outputs `outcvA` and `outcvb` are modified in-place to avoid memory
reallocation."""
# GELEES_riemann_solver_setup( self.mesh.C, self.mesh.nC, self.W, self.cvA, self.cvb)
cvC: FloatArray = self.mesh.C
cvnC: FloatArray = self.mesh.nC
#cvZ: FloatArray = self.cvZ
csig: FloatArray = self.get_primal_variable("sigma")
cu : FloatArray = self.get_primal_variable("u")
crho = self.get_primal_variable("rho")
lamb = self.law.lamb
mu = self.law.mu
for j in range(self.cvA.shape[0]):
rhoj = crho[j]
#print(rhoj)
aLj = np.sqrt((lamb+2*mu)/rhoj)
aTj = np.sqrt(mu/rhoj)
vCj = cvC[j]
vnCj = cvnC[j]
for r in range(self.cvA.shape[1]):
Cjr = vCj[r]
nCjr = vnCj[r]
CjrxCjr_00 = Cjr[0]*Cjr[0]
CjrxCjr_01 = Cjr[0]*Cjr[1]
CjrxCjr_10 = CjrxCjr_01
CjrxCjr_11 = Cjr[1]*Cjr[1]
self.cvA[j,r][0,0] = rhoj *aLj * CjrxCjr_00/nCjr +rhoj*aTj*( nCjr - CjrxCjr_00/nCjr )
self.cvA[j,r][0,1] = rhoj *aLj * CjrxCjr_01/nCjr - rhoj*aTj*( CjrxCjr_01/nCjr )
self.cvA[j,r][1,0] = rhoj *aLj * CjrxCjr_10/nCjr - rhoj*aTj*( CjrxCjr_10/nCjr )
self.cvA[j,r][1,1] = rhoj *aLj * CjrxCjr_11/nCjr +rhoj*aTj*( nCjr - CjrxCjr_11/nCjr )
#self.cvA[:]= cvZ*cvnC[..., np.newaxis, np.newaxis]
# b_jr = sigj . Cjr + A_jr . u
csig = csig[:,np.newaxis,...]
cvC = cvC[..., np.newaxis]
self.cvb[:] = -(csig@cvC).squeeze()
self.cvb[:] += self.cvA[...,0]*cu[:, np.newaxis, np.newaxis, 0] \
+ self.cvA[...,1]*cu[:, np.newaxis, np.newaxis, 1]
def _riemann_solver_setup_1(self) -> None :
lamb = self.law.lamb
mu = self.law.mu
cvC: FloatArray = self.mesh.C
cvnC: FloatArray = self.mesh.nC
crho, csig, cu = self.get_primal_variable(["rho", "sigma", "u"])
cvA, cvb = self.cvA, self.cvb
riemann_setup_1_guvec(lamb, mu, crho, cvC, cvnC, csig, cu, cvA, cvb)
def _riemann_solver_setup(self)-> None:
r"""
Computes the local cell-vertex matrix and vector contributions for the
nodal acoustic Rie
mann solver.
This function vectorizes the calculation of the $A_{jr}$ matrices and
$b_{jr}$ vectors for all cells and their vertices simultaneously. These
local contributions will later be assembled into a global linear system
at each node to compute the nodal velocities.
The continuous formulations for the acoustic solver are:
$A_{jr} = \rho_j c_j \frac{C_{jr} \otimes C_{jr}}{|C_{jr}|}$
$b_{jr} = p_j C_{jr} + A_{jr} \cdot u_j$
Parameters
----------
cvC : FloatArray
3D array of shape `(ncells, 3, 2)` containing the corner normal vectors
$C_{jr}$ for each vertex of each cell.
cvnC : FloatArray
2D array of shape `(ncells, 3)` containing the Euclidean norms (magnitudes)
of the corner normal vectors $|C_{jr}|$.
W : FloatArray
2D array of shape `(ncells, 5)` containing the primitive variables for
each cell (density, x-velocity, y-velocity, pressure, sound speed).
outcvA : FloatArray
4D output array of shape `(ncells, 3, 2, 2)`. The computed symmetric
$2 \times 2$ matrices $A_{jr}$ are stored here. Modified in-place.
outcvb : FloatArray
3D output array of shape `(ncells, 3, 2)`. The computed 2D vectors
$b_{jr}$ are stored here. Modified in-place.
Returns
-------
None
The outputs `outcvA` and `outcvb` are modified in-place to avoid memory
reallocation."""
# GELEES_riemann_solver_setup( self.mesh.C, self.mesh.nC, self.W, self.cvA, self.cvb)
cvC: FloatArray = self.mesh.C
cvnC: FloatArray = self.mesh.nC
cvZ: FloatArray = self.cvZ
cvA: FloatArray = self.cvA
cvb: FloatArray = self.cvb
csig: FloatArray = self.get_primal_variable("sigma")
cu : FloatArray = self.get_primal_variable("u")
if USE_NUMBA:
GlaceExplicitLagrangian_HyperElasticty_EquationSolver._riemann_solver_setup_njit(cu, csig, cvZ, cvC, cvnC, cvA, cvb)
else:
cvA[:]= cvZ*cvnC[..., np.newaxis, np.newaxis]
# b_jr = sigj . Cjr + A_jr . u
csig = csig[:,np.newaxis,...]
cvC = cvC[..., np.newaxis]
cvb[:] = -(csig@cvC).squeeze()
cvb[:] += 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, csig, cvZ, cvC, cvnC, out_cvA, out_cvb):
N =cu.shape[0]
for j in numba.prange(N):
sigj = csig[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] = -sigj[0,0]* Cjr[0] -sigj[0,1]*Cjr[1] + Ajr[0, 0]*uj0 + Ajr[0, 1]*uj1
bjr[1] = -sigj[1,0]* Cjr[0] -sigj[1,1]*Cjr[1] + Ajr[1, 0]*uj0 + Ajr[1, 1]*uj1
[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).
"""
cV = self.mesh.cell_meas
#cc = self.get_primal_variable("c") # per cell sound speed
cc = self.cc
cnC = self.mesh.nC
return type(self).critical_timestep_pnjit(cV, cnC, cc)
[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
def update_conservative():
update_conservative_vars_pnjit(cvqflux, cvu, dt, self.conservation_law.U)
update_conservative()
def update_xy():
#mesh.xy[:] += dt*vu
axpy(dt, vu, mesh.xy, parallel =False )
update_xy()
# update primal variable (on moved cells)
update_MeshStruct(mesh.xy, mesh)
update_AdditionalMeshStruct(mesh, self.mesh_adata)
#m = self.conservation_law.U[:, I_M]
m, mF = self.get_conservative_variable(["m", "mF"])
np.multiply(self.mesh_adata.F, m[..., None, None], out = mF)
self.conservation_law.conservative_to_primal( mesh.cell_meas)
self.compute_impedance_tensors(self.mesh, self.cvZ)
GlaceExplicitLagrangian_HyperElasticty_EquationSolver.cell_soundspeed(self.cvZ, self.cc)
#update time
self.t += dt
#compute critical time step
dtcnew = self.critical_timestep()
return dtcnew
step = step_np