Merge branch 'ec-34287_DYNA_ISS_VARI_opti' into 'main' (60f41f0d) · Commits · codeaster / src

code_aster/MacroCommands/Miss/calc_miss_vari.py

+81 −49

Original line number	Diff line number	Diff line
		@@ -22,12 +22,13 @@ from math import pi, sqrt
		import aster
		import numpy as NP
		from numpy import linalg
		from typing import List

		from ...Cata.Syntax import _F
		from ...CodeCommands import COMB_MATR_ASSE, CREA_CHAMP, DYNA_VIBRA, LIRE_FORC_MISS, LIRE_IMPE_MISS
		from ...Messages import UTMESS
		from ...Utilities import disable_fpe
		from ..Utils.signal_correlation_utils import CALC_COHE
		from ..Utils.signal_correlation_utils import calc_coherency_matrix


		def calc_miss_vari(self):
		@@ -42,17 +43,23 @@ def calc_miss_vari(self):
		list_NOM_CMP = [self.NOM_CMP]
		self.NOM_CMP = None

		# BOUCLE SUR LES FREQUENCES : Matrice de cohérence
		PVEC = [None] * NB_FREQ
		nbm = [None] * NB_FREQ
		for k in range(NB_FREQ):
		freqk = self.FREQ_INIT + self.FREQ_PAS * k
		if self.INFO == 2:
		aster.affiche("MESSAGE", "FREQUENCE DE CALCUL: " + str(freqk))
		# COHERENCE
		COHE = CALC_COHE(freqk * 2.0 * pi, **self.cohe_params)
		# EIGENVALUE DECOMP
		PVEC[k], nbm[k] = compute_POD(COHE, self.calc_params["PRECISION"], self.INFO)
		frequencies = self.FREQ_INIT + self.FREQ_PAS * NP.arange(NB_FREQ)
		kwargs = self.cohe_params

		coherency_matrix = calc_coherency_matrix(
		frequencies=frequencies,
		model=kwargs["TYPE"],
		nom_mail=kwargs["MAILLAGE"],
		nom_group_inter=kwargs["GROUP_NO_INTERF"],
		**kwargs,
		)

		PVEC = compute_pod(
		coherency_matrix=coherency_matrix,
		precision=self.calc_params["PRECISION"],
		frequencies=frequencies,
		info=self.INFO,
		)

		# RECUPERATION DES MODES MECA (STATIQUES)
		nbmodd = self.mat_gene_params["NBMODD"]
		@@ -69,7 +76,7 @@ def calc_miss_vari(self):
		dict_modes["NOM_CMP"] = nom_cmp
		# BOUCLE SUR LES FREQUENCES
		for k in range(NB_FREQ):
		dict_modes["nbpod"] = nbm[k]
		dict_modes["nbpod"] = len(PVEC[k])
		# CALCUL ISS VARI
		if self.interf_params["MODE_INTERF"] != "QUELCONQUE":
		RESU[i_cmp] = self.compute_freqk(k, RESU[i_cmp], PVEC[k], dict_modes)
		@@ -85,42 +92,67 @@ def calc_miss_vari(self):
		# --------------------------------------------------------------------------------
		# ROUTINES
		# --------------------------------------------------------------------------------
		def compute_POD(COHE, PRECISION, INFO):
		"""compute POD"""
		# EIGENVALUE DECOMP
		# nbno = self.interf_params['NBNO']
		# On desactive temporairement les FPE
		def compute_pod(
		coherency_matrix: NP.ndarray, precision: float, frequencies: NP.ndarray, info: int
		) -> List[NP.ndarray]:
		"""compute POD
		Args:
		coherency_matrix : coherency matrix shape (nb_freq, nb_nodes, nb_nodes)
		precision : unitary threshold of the sum of the energy
		frequencies : frequencies
		info : verbosity level

		"""

		def iter_on_eigen_values_and_vectors_all():
		"""Resolve once the the eigenvalue problem for all frequencies"""
		if info == 2:
		aster.affiche("MESSAGE", "RESOLUTION DES VALEURS PROPRES POUR TOUTES LES FREQUENCES")

		with disable_fpe():
		eig, vec = linalg.eig(COHE)
		vec = NP.transpose(vec) # les vecteurs sont en colonne dans numpy
		eig, vec = linalg.eigh(coherency_matrix)

		eig = eig.real
		vec = vec.real
		# on rearrange selon un ordre decroissant
		eig = NP.where(eig < 1.0e-10, 0.0, eig)
		order = NP.argsort(eig)[::-1]
		eig = NP.take(eig, order)
		vec = NP.take(vec, order, 0)
		# Nombre de modes POD a retenir
		etot = NP.sum(eig**2)
		ener = 0.0
		nbme = 0
		nbno = len(eig)
		while nbme < nbno:
		ener = eig[nbme] ** 2 + ener
		prec = ener / etot
		nbme = nbme + 1
		if INFO == 2:
		texte = "VALEUR PROPRE " + str(nbme) + " : " + str(eig[nbme - 1])
		aster.affiche("MESSAGE", texte)
		if prec > PRECISION:

		eig[eig < 1.0e-10] = 0.0
		# les vecteurs sont en colonne en sortie de linalg.eigh
		vec = NP.transpose(vec, axes=(0, 2, 1))
		# Sort eigenvalues and eigenvectors in descending order
		eig = eig[:, ::-1] # shape (nb_fre, nb_nodes)
		vec = vec[:, ::-1]
		yield from zip(eig, vec)

		pvec = []

		for freq, (values, vectors) in zip(frequencies, iter_on_eigen_values_and_vectors_all()):
		cum_nrj = 0.0

		energies = values**2
		e_tot = NP.sum(energies)

		if info == 2:
		aster.affiche("MESSAGE", "FREQUENCE DE CALCUL: " + str(freq))

		# keep only the modes than explain the precision of the cumulative energy
		for nbme, (nrj, value) in enumerate(zip(energies, values), 1):
		if info == 2:
		aster.affiche("MESSAGE", "VALEUR PROPRE " + str(nbme) + " : " + str(value))

		cum_nrj += nrj
		prec = cum_nrj / e_tot

		if prec > precision:
		break
		if INFO == 2:

		if info == 2:
		aster.affiche("MESSAGE", "NOMBRE DE MODES POD RETENUS : " + str(nbme))
		aster.affiche("MESSAGE", "PRECISION (ENERGIE RETENUE) : " + str(prec))
		PVEC = NP.zeros((nbme, nbno))
		for k1 in range(nbme):
		PVEC[k1, 0:nbno] = NP.sqrt(eig[k1]) * vec[k1]
		return PVEC, nbme

		_pvec = NP.sqrt(values[:nbme, None]) * vectors[:nbme]
		pvec.append(_pvec)

		return pvec


		# ---------------------------------------------------------------------

code_aster/MacroCommands/Utils/calc_coherency.py

+1 −77

Original line number	Diff line number	Diff line
		# coding=utf-8
		# --------------------------------------------------------------------
		# Copyright (C) 1991 - 2022 - EDF R&D - www.code-aster.org
		# Copyright (C) 1991 - 2025 - EDF R&D - www.code-aster.org
		# This file is part of code_aster.
		#
		# code_aster is free software: you can redistribute it and/or modify
		@@ -95,79 +95,3 @@ def calc_cohefromdata(acce1, acce2, dt, Mm):
		gami2 = Sxf12[i2] / np.sqrt(np.array(Sxf1[i2]) * np.array(Sxf2[i2]))
		gamw.append(gami2)
		return wm / 2.0 / pi, (np.mean(gamw, 0))


		# -------------------------------------------------------------------
		# COHERENCY MATRIX
		# --------------------------------------------------------------------
		def CALC_COHE(freq, **kwargs):
		# Frequency is in rad/s: freq= f2pi
		# kwargs: VITE_ONDE, PARA_ALPHA, TYPE, MAILLAGE,
		model = kwargs["TYPE"]
		nom_mail = kwargs["MAILLAGE"]
		nom_group_inter = kwargs["GROUP_NO_INTERF"]
		if "NOEUDS_INTERF" in kwargs:
		noe_interf = kwargs["NOEUDS_INTERF"]
		else:
		liste_nom, noe_interf = get_group_nom_coord(nom_group_inter, nom_mail)
		if "DIST" in kwargs:
		DIST2 = kwargs["DIST"]
		else:
		DIST2 = calc_dist2(noe_interf)
		# # ----MITA & LUCO
		if model == "MITA_LUCO":
		# PARAMETRES fonction de coherence
		VITE_ONDE = kwargs["VITE_ONDE"]
		alpha = kwargs["PARA_ALPHA"]
		COHE = NP.exp(-(DIST2 * (alpha * freq / VITE_ONDE) ** 2.0))
		# ----ABRAHAMSON ROCK (EPRI)
		elif model == "ABRAHAMSON":
		p_a1 = 1.647
		p_a2 = 1.01
		p_a3 = 0.4
		p_n1 = 7.02
		nbno = len(noe_interf)
		freqk = freq / (2.0 * pi)
		COHE = NP.zeros((nbno, nbno))
		for no1 in range(nbno):
		for no2 in range(nbno):
		# dist_xi = sqrt((XX[no1] - XX[no2])2 + (YY[no1] - YY[no2])2)
		dist_xi = sqrt(DIST2[no1, no2])
		p_n2 = 5.1 - 0.51 * log(dist_xi + 10.0)
		pfc = -1.886 + 2.221 * log(4000.0 / (dist_xi + 1.0) + 1.5)
		term1 = 1.0 + (freqk * tanh(p_a3 * dist_xi) / (p_a1 * pfc)) ** p_n1
		term2 = 1.0 + (freqk * tanh(p_a3 * dist_xi) / (p_a2 * pfc)) ** p_n2
		COHE[no1, no2] = 1.0 / sqrt(term1 * term2)
		elif model == "ABRA_ROCHER":
		p_a1 = 1.0
		p_a2 = 40.0
		p_a3 = 0.4
		p_n2 = 16.4
		nbno = len(noe_interf)
		freqk = freq / (2.0 * pi)
		COHE = NP.zeros((nbno, nbno))
		for no1 in range(nbno):
		for no2 in range(nbno):
		dist_xi = sqrt(DIST2[no1, no2])
		p_n1 = 3.8 - 0.04 * log(dist_xi + 1.0) + 0.0105 * (log(dist_xi + 1.0) - 3.6) ** 2
		pfc = 27.9 - 4.82 * log(dist_xi + 1.0) + 1.24 * (log(dist_xi + 1.0) - 3.6) ** 2
		term1 = 1.0 + (freqk * tanh(p_a3 * dist_xi) / (p_a1 * pfc)) ** p_n1
		term2 = 1.0 + (freqk * tanh(p_a3 * dist_xi) / (p_a2)) ** p_n2
		COHE[no1, no2] = 1.0 / sqrt(term1 * term2)
		elif model == "ABRA_SOLMOYEN":
		p_a1 = 1.0
		p_a3 = 0.4
		p_n1 = 3.0
		p_n2 = 15.0
		nbno = len(noe_interf)
		freqk = freq / (2.0 * pi)
		COHE = NP.zeros((nbno, nbno))
		for no1 in range(nbno):
		for no2 in range(nbno):
		dist_xi = sqrt(DIST2[no1, no2])
		p_a2 = 15.8 - 0.044 * dist_xi
		pfc = 14.3 - 2.35 * log(dist_xi + 1.0)
		term1 = 1.0 + (freqk * tanh(p_a3 * dist_xi) / (p_a1 * pfc)) ** p_n1
		term2 = 1.0 + (freqk * tanh(p_a3 * dist_xi) / (p_a2)) ** p_n2
		COHE[no1, no2] = 1.0 / sqrt(term1 * term2)
		return COHE

code_aster/MacroCommands/Utils/signal_correlation_utils.py

+55 −43

Original line number	Diff line number	Diff line
		@@ -47,12 +47,17 @@ from .random_signal_utils import ACCE2SROM, acce_filtre_CP, calc_dsp_FR, calc_ds
		# -------------------------------------------------------------------
		# COHERENCY MATRIX
		# --------------------------------------------------------------------
		def CALC_COHE(freq, **kwargs):
		# Frequency is in rad/s: freq= f2pi
		# kwargs: VITE_ONDE, PARA_ALPHA, TYPE, MAILLAGE,
		model = kwargs["TYPE"]
		nom_mail = kwargs["MAILLAGE"]
		nom_group_inter = kwargs["GROUP_NO_INTERF"]
		def calc_coherency_matrix(frequencies, model, nom_mail, nom_group_inter, **kwargs):
		"""
		Calculation of coherency matrix

		Arguments:
		frequencies : list of frequencies
		model : coherency function type
		nom_mail : mesh name
		nom_group_inter : group of nodes constituting the soil-structure interface
		kwargs: other arguments from the commands
		"""
		if "NOEUDS_INTERF" in kwargs:
		noe_interf = kwargs["NOEUDS_INTERF"]
		else:
		@@ -61,65 +66,72 @@ def CALC_COHE(freq, **kwargs):
		DIST2 = kwargs["DIST"]
		else:
		DIST2 = calc_dist2(noe_interf)
		# # ----MITA & LUCO
		frequencies = frequencies[:, None, None]
		dist_xi = NP.sqrt(DIST2)
		if model == "MITA_LUCO":
		# PARAMETRES fonction de coherence
		VITE_ONDE = kwargs["VITE_ONDE"]
		alpha = kwargs["PARA_ALPHA"]
		COHE = NP.exp(-(DIST2 * (alpha * freq / VITE_ONDE) ** 2.0))
		# ----ABRAHAMSON ROCK (EPRI)
		COHE = NP.exp(-(DIST2 * (alpha * frequencies * 2.0 * pi / VITE_ONDE) ** 2.0))
		elif model == "ABRAHAMSON":
		p_a1 = 1.647
		p_a2 = 1.01
		p_a3 = 0.4
		p_n1 = 7.02
		nbno = len(noe_interf)
		freqk = freq / (2.0 * pi)
		COHE = NP.zeros((nbno, nbno))
		for no1 in range(nbno):
		for no2 in range(nbno):
		# dist_xi = sqrt((XX[no1] - XX[no2])2 + (YY[no1] - YY[no2])2)
		dist_xi = sqrt(DIST2[no1, no2])
		p_n2 = 5.1 - 0.51 * log(dist_xi + 10.0)
		pfc = -1.886 + 2.221 * log(4000.0 / (dist_xi + 1.0) + 1.5)
		term1 = 1.0 + (freqk * tanh(p_a3 * dist_xi) / (p_a1 * pfc)) ** p_n1
		term2 = 1.0 + (freqk * tanh(p_a3 * dist_xi) / (p_a2 * pfc)) ** p_n2
		COHE[no1, no2] = 1.0 / sqrt(term1 * term2)
		p_n2 = 5.1 - 0.51 * NP.log(dist_xi + 10.0)
		pfc = -1.886 + 2.221 * NP.log(4000.0 / (dist_xi + 1.0) + 1.5)
		numerator = frequencies * NP.tanh(p_a3 * dist_xi)
		term1 = 1.0 + (numerator / (p_a1 * pfc)) ** p_n1
		term2 = 1.0 + (numerator / (p_a2 * pfc)) ** p_n2
		COHE = 1.0 / NP.sqrt(term1 * term2)
		elif model == "ABRA_ROCHER":
		p_a1 = 1.0
		p_a2 = 40.0
		p_a3 = 0.4
		p_n2 = 16.4
		nbno = len(noe_interf)
		freqk = freq / (2.0 * pi)
		COHE = NP.zeros((nbno, nbno))
		for no1 in range(nbno):
		for no2 in range(nbno):
		dist_xi = sqrt(DIST2[no1, no2])
		p_n1 = 3.8 - 0.04 * log(dist_xi + 1.0) + 0.0105 * (log(dist_xi + 1.0) - 3.6) ** 2
		pfc = 27.9 - 4.82 * log(dist_xi + 1.0) + 1.24 * (log(dist_xi + 1.0) - 3.6) ** 2
		term1 = 1.0 + (freqk * tanh(p_a3 * dist_xi) / (p_a1 * pfc)) ** p_n1
		term2 = 1.0 + (freqk * tanh(p_a3 * dist_xi) / (p_a2)) ** p_n2
		COHE[no1, no2] = 1.0 / sqrt(term1 * term2)
		log_d = NP.log(dist_xi + 1.0)
		p_n1 = 3.8 - 0.04 * log_d + 0.0105 * (log_d - 3.6) ** 2
		pfc = 27.9 - 4.82 * log_d + 1.24 * (log_d - 3.6) ** 2
		numerator = frequencies * NP.tanh(p_a3 * dist_xi)
		term1 = 1.0 + (numerator / (p_a1 * pfc)) ** p_n1
		term2 = 1.0 + (numerator / p_a2) ** p_n2
		COHE = 1.0 / NP.sqrt(term1 * term2)
		elif model == "ABRA_SOLMOYEN":
		p_a1 = 1.0
		p_a3 = 0.4
		p_n1 = 3.0
		p_n2 = 15.0
		nbno = len(noe_interf)
		freqk = freq / (2.0 * pi)
		COHE = NP.zeros((nbno, nbno))
		for no1 in range(nbno):
		for no2 in range(nbno):
		dist_xi = sqrt(DIST2[no1, no2])
		p_a2 = 15.8 - 0.044 * dist_xi
		pfc = 14.3 - 2.35 * log(dist_xi + 1.0)
		term1 = 1.0 + (freqk * tanh(p_a3 * dist_xi) / (p_a1 * pfc)) ** p_n1
		term2 = 1.0 + (freqk * tanh(p_a3 * dist_xi) / (p_a2)) ** p_n2
		COHE[no1, no2] = 1.0 / sqrt(term1 * term2)
		pfc = 14.3 - 2.35 * NP.log(dist_xi + 1.0)

		numerator = frequencies * NP.tanh(p_a3 * dist_xi)
		term1 = 1.0 + (numerator / (p_a1 * pfc)) ** p_n1
		term2 = 1.0 + (numerator / p_a2) ** p_n2
		COHE = 1.0 / NP.sqrt(term1 * term2)
		return COHE


		def CALC_COHE(puls, **kwargs):
		"""
		Use to call calc_coherency_matrix from CALC_MISS
		Convert pulsations to frequencies and call calc_coherency_matrix

		Arguments :
		pusl : pulsation values vector

		Return :
		Coherency matrix
		"""
		coherency_matrix = calc_coherency_matrix(
		frequencies=NP.array([puls / (2.0 * pi)]),
		model=kwargs["TYPE"],
		nom_mail=kwargs["MAILLAGE"],
		nom_group_inter=kwargs["GROUP_NO_INTERF"],
		**kwargs
		)
		return coherency_matrix[0]


		def get_group_nom_coord(group_inter, nom_mail):
		print("in signal_correlation_utils")
		# no des noeuds