MMGP_demo / utils_inference.py
fabiencasenave's picture
initial commit
f4d7da3
raw
history blame contribute delete
No virus
12.5 kB
import Muscat.Containers.ElementsDescription as ED
from Muscat.FE.FETools import PrepareFEComputation
from Muscat.FE.Fields.FEField import FEField
import Muscat.Containers.MeshModificationTools as MMT
from Muscat.Containers import MeshGraphTools as MGT
from Muscat.Containers.Filters import FilterObjects as FO
from Muscat.Containers.Filters import FilterOperators as FOp
import Muscat.Containers.MeshInspectionTools as MIT
out_fields_names = ['Ux', 'Uy', 'p', 'nut']
in_scalars_names = ['angle_of_attack', 'inlet_velocity']
from Muscat.Containers.NativeTransfer import NativeTransfer
import copy
import bisect
import numpy as np
from scipy import spatial
import pickle
from Muscat.Bridges.CGNSBridge import CGNSToMesh
import torch
device = torch.device("cpu")
num_steps = [2_000]
from plaid.containers.dataset import Sample
def pretreat_sample(mesh):
###################
# Compute the skin of the mesh (containing the external boundary and the airfoil boundary)
###################
MMT.ComputeSkin(mesh, md = 2, inPlace = True)
###################
# Extract ids of the bar elements corresponding to the airfoil
###################
ff1 = FO.ElementFilter(zone = lambda p: (-p[:,0]-1.99))
ff2 = FO.ElementFilter(zone = lambda p: (p[:,0]-3.99))
ff3 = FO.ElementFilter(zone = lambda p: (-p[:,1]-1.49))
ff4 = FO.ElementFilter(zone = lambda p: (p[:,1]-1.49))
efAirfoil = FOp.IntersectionFilter(filters=[ff1, ff2, ff3, ff4])
airfoil_ids = efAirfoil.GetIdsToTreat(mesh, "bar2")
###################
# Preparations
###################
ext_bound = np.setdiff1d(mesh.elements["bar2"].GetTag("Skin").GetIds(), airfoil_ids)
mesh.elements["bar2"].GetTag("External_boundary").SetIds(ext_bound)
mesh.elements["bar2"].GetTag("Airfoil").SetIds(airfoil_ids)
nfExtBoundary = FO.NodeFilter(eTag = "External_boundary")
nodeIndexExtBoundary = nfExtBoundary.GetNodesIndices(mesh)
mesh.GetNodalTag("External_boundary").AddToTag(nodeIndexExtBoundary)
nfAirfoil = FO.NodeFilter(eTag = "Airfoil")
nodeIndexAirfoil = nfAirfoil.GetNodesIndices(mesh)
mesh.GetNodalTag("Airfoil").AddToTag(nodeIndexAirfoil)
###################
# Add node tag for the intrado and extrado
###################
nfAirfoil = FO.NodeFilter(eTag = "Airfoil")
nodeIndexAirfoil = nfAirfoil.GetNodesIndices(mesh)
mesh.GetNodalTag("Airfoil").AddToTag(nodeIndexAirfoil)
airfoil = ExtractAirfoil(mesh)
indices_extrado = airfoil[0][0]
indices_intrado = airfoil[0][1]
mesh.GetNodalTag("Extrado").AddToTag(indices_extrado)
mesh.GetNodalTag("Intrado").AddToTag(indices_intrado)
efExtrado = FO.ElementFilter(nTag = "Extrado")
efIntrado = FO.ElementFilter(nTag = "Intrado")
mesh.elements["bar2"].GetTag("Extrado").SetIds(efExtrado.GetIdsToTreat(mesh, "bar2"))
mesh.elements["bar2"].GetTag("Intrado").SetIds(efIntrado.GetIdsToTreat(mesh, "bar2"))
def ExtractPathFromMeshOfBars(mesh, startingClosestToPoint, trigo_dir = True):
nodeGraph0Airfoild = MGT.ComputeNodeToNodeGraph(mesh, dimensionality=1)
nodeGraphAirfoild = [list(nodeGraph0Airfoild[i].keys()) for i in range(nodeGraph0Airfoild.number_of_nodes())]
tree = spatial.KDTree(mesh.nodes)
_, indicesTrailEdge = tree.query([startingClosestToPoint], k=1)
p1init = indicesTrailEdge[0]
temp1=mesh.nodes[nodeGraphAirfoild[p1init][0]][1]
temp2=mesh.nodes[nodeGraphAirfoild[p1init][1]][1]
if trigo_dir:
condition = temp1 > temp2
else:
condition = temp1 < temp2
if condition:
p2 = nodeGraphAirfoild[p1init][0]
else:
p2 = nodeGraphAirfoild[p1init][1]
p1 = p1init
path = [p1, p2]
while p2 != p1init:
p2save = p2
tempArray = np.asarray(nodeGraphAirfoild[p2])
p2 = tempArray[tempArray!=p1][0]
p1 = p2save
path.append(p2)
return path
def ExtractAirfoil(mesh):
efAirfoil = FO.ElementFilter(elementType=ED.Bar_2, eTag=["Airfoil"])
airfoilMesh = MIT.ExtractElementsByElementFilter(mesh, efAirfoil)
path = ExtractPathFromMeshOfBars(airfoilMesh, np.array([1.,0.]))
tree = spatial.KDTree(airfoilMesh.nodes[path])
_, indicesLeadEdge = tree.query([[0.,0.]], k=1)
indices_extrado = path[:indicesLeadEdge[0]+1]
indices_intrado = path[indicesLeadEdge[0]:]
indices_airfoil = [indices_extrado, indices_intrado]
nodes_extrado = mesh.nodes[indices_extrado]
nodes_intrado = mesh.nodes[indices_intrado]
nodes_airfoil = [nodes_extrado, nodes_intrado]
return indices_airfoil, nodes_airfoil
def computeAirfoilCurvAbscissa(airfoil):
indices_airfoil = airfoil[0]
nodes_airfoil = airfoil[1]
curv_abscissa = []
for i in range(2):
local_curv_abscissa = np.zeros(len(indices_airfoil[i]))
for j in range(1,len(local_curv_abscissa)):
local_curv_abscissa[j] = local_curv_abscissa[j-1] + np.linalg.norm(nodes_airfoil[i][j]-nodes_airfoil[i][j-1])
local_curv_abscissa /= local_curv_abscissa[-1]
curv_abscissa.append(local_curv_abscissa)
return curv_abscissa
def MapAirfoil(airfoil_ref, curv_abscissa_ref, curv_abscissa):
nodes_airfoil_ref = airfoil_ref[1]
dim_nodes = nodes_airfoil_ref[0][0].shape[0]
mapped_airfoil = []
for i in range(2):
local_mapped_airfoil = np.zeros((len(curv_abscissa[i])-1, dim_nodes))
for j in range(len(curv_abscissa[i])-1):
index = max(bisect.bisect_right(curv_abscissa_ref[i], curv_abscissa[i][j]) - 1, 0)
a = nodes_airfoil_ref[i][index]
b = nodes_airfoil_ref[i][index+1]
dl = curv_abscissa[i][j] - curv_abscissa_ref[i][index]
dir = (b-a)/np.linalg.norm(b-a)
local_mapped_airfoil[j] = a + dl * dir
mapped_airfoil.append(local_mapped_airfoil)
return mapped_airfoil
def GetFieldTransferOpCppStep1(inputField, nbThreads = None):
method="Interp/Clamp"
nt = NativeTransfer()
if nbThreads is not None:
nt.SetMaxConcurrency(nbThreads)
nt.SetTransferMethod(method)
defaultOptions = {"usePointSearch": True,
"useElementSearch": False,
"useElementSearchFast": False,
"useEdgeSearch": True,
}
options = {}
defaultOptions.update(options)
dispatch = {"usePointSearch": nt.SetUsePointSearch,
"useElementSearch": nt.SetUseElementSearch,
"useElementSearchFast": nt.SetUseElementSearchFast,
"useEdgeSearch": nt.SetUseEdgeSearch,
"DifferentialOperator": nt.SetDifferentialOperator,
}
for k, v in defaultOptions.items():
if k in dispatch.keys():
dispatch[k](v)
else:
raise RuntimeError(f"Option {k} not valid")
nt.SetSourceFEField(inputField, None)
return nt
def GetFieldTransferOpCppStep2(nt, targetPoints):
nt.SetTargetPoints(targetPoints)
nt.Compute()
op = nt.GetOperator()
status = nt.GetStatus()
return op, status
def morph_sample(mesh, airfoil_0, curv_abscissa_0):
airfoil_1 = ExtractAirfoil(mesh)
curv_abscissa_1 = computeAirfoilCurvAbscissa(airfoil_1)
mapped_airfoil = MapAirfoil(airfoil_0, curv_abscissa_0, curv_abscissa_1)
##############################################################
# Compute global target displacement and masks for RBF field morphing
##############################################################
indices_extrado_to_morph_1 = airfoil_1[0][0][:-1]
indices_intrado_to_morph_1 = airfoil_1[0][1][:-1]
other_boundary_ids_1 = mesh.GetNodalTag("External_boundary").GetIds()
l1 = len(indices_extrado_to_morph_1)
l2 = len(indices_intrado_to_morph_1)
l3 = len(other_boundary_ids_1)
targetDisplacement = np.zeros((l1 + l2 + l3, 2))
targetDisplacementMask = np.zeros((l1 + l2 + l3), dtype = int)
targetDisplacement[:l1] = mapped_airfoil[0] - mesh.nodes[indices_extrado_to_morph_1]
targetDisplacement[l1:l1+l2] = mapped_airfoil[1] - mesh.nodes[indices_intrado_to_morph_1]
targetDisplacement[l1+l2:l1+l2+l3] = np.zeros((l3,2))
targetDisplacementMask[:l1] = indices_extrado_to_morph_1
targetDisplacementMask[l1:l1+l2] = indices_intrado_to_morph_1
targetDisplacementMask[l1+l2:l1+l2+l3] = other_boundary_ids_1
##############################################################
# Compute the morphing
##############################################################
# RBF morphing
mesh_nodes = mesh.nodes.copy()
morphed_nodes = MMT.Morphing(mesh, targetDisplacement, targetDisplacementMask, radius=None)
mesh.nodes = morphed_nodes
mesh.nodeFields['X'] = mesh_nodes[:,0]
mesh.nodeFields['Y'] = mesh_nodes[:,1]
return mesh
def project_sample(morphed_mesh, morphed_mesh_0):
projected_mesh = copy.deepcopy(morphed_mesh_0)
space_, numberings_, _, _ = PrepareFEComputation(morphed_mesh)
inputFEField = FEField(name="dummy", mesh=morphed_mesh, space=space_, numbering=numberings_[0])
nt = GetFieldTransferOpCppStep1(inputFEField, 1)
FE_interpolation_op, _ = GetFieldTransferOpCppStep2(nt, morphed_mesh_0.nodes)
for pfn in out_fields_names + ['X', 'Y']:
projected_mesh.nodeFields[pfn] = FE_interpolation_op.dot(morphed_mesh.nodeFields[pfn])
return projected_mesh
def pretreat_morph_and_project_mesh(i_sample, dataset, indices, morphed_mesh_0, airfoil_0, curv_abscissa_0):
###################
## 1) Pretreat data
###################
sample = Sample.model_validate(pickle.loads(dataset[int(indices[i_sample])]["sample"]))
mesh = CGNSToMesh(sample.get_mesh())
pretreat_sample(mesh)
###################
## 2) Morph data
###################
morphed_mesh = morph_sample(mesh, airfoil_0, curv_abscissa_0)
# print("morphed_mesh =", morphed_mesh)
# from Muscat.IO import XdmfWriter as XW
# XW.WriteMeshToXdmf("morphed_mesh.xdmf", morphed_mesh)
###################
## 3) Project data
###################
projected_mesh = project_sample(morphed_mesh, morphed_mesh_0)
# print("projected_mesh =", projected_mesh)
# from Muscat.IO import XdmfWriter as XW
# XW.WriteMeshToXdmf("projected_mesh.xdmf", projected_mesh)
in_scalars = [sample.get_scalar(sn) for sn in in_scalars_names]
return [projected_mesh, in_scalars, morphed_mesh.nodes]
def infer(dataset, indice, training_data):
common_mesh, correlationOperator2c, airfoil_0, curv_abscissa_0, scalar_scalers, scalerX, pca_clouds, pca_fields, y_scalers, kmodels, X_train, y_train = training_data
projected_mesh, in_scalars, morphed_mesh_nodes = pretreat_morph_and_project_mesh(indice, dataset = dataset, indices = np.arange(len(dataset)), morphed_mesh_0 = common_mesh, airfoil_0 = airfoil_0, curv_abscissa_0 = curv_abscissa_0)
space_, numberings_, _, _ = PrepareFEComputation(common_mesh)
inputFEField_0 = FEField(name="dummy", mesh=common_mesh, space=space_, numbering=numberings_[0])
nt_0 = GetFieldTransferOpCppStep1(inputFEField_0, nbThreads = 2)
iffo = GetFieldTransferOpCppStep2(nt_0, morphed_mesh_nodes)[0]
clouds = np.stack([projected_mesh.nodeFields["X"], projected_mesh.nodeFields["Y"]], axis=1)
scalars = np.array(in_scalars)
clouds = clouds.reshape(-1, 1)
X_pca = np.dot(pca_clouds, correlationOperator2c.dot(clouds)).T
X_scalars = scalar_scalers.transform(scalars.reshape(1, -1))
unscaled_X = np.concatenate([X_pca, X_scalars], axis=-1)
X = scalerX.transform(unscaled_X)
predictions = {}
y = []
for i, fn in enumerate(out_fields_names):
X_ = torch.tensor(X, dtype=torch.float64).to(device)
output_dim = pca_fields[0].shape[0]
n_samples = 1
y_pred_i = np.empty((n_samples, output_dim, len(num_steps)))
for j in range(output_dim):
for k in range(len(num_steps)):
y_pred_i[:,j,k] = kmodels[i][j][k](X_).detach().cpu().numpy().squeeze()
y.append(y_pred_i)
y_pred_i = np.mean(y_pred_i, axis = 2)
y_pred_i = y_scalers[i].inverse_transform(y_pred_i)
y_pred_common_i = np.dot(y_pred_i, pca_fields[i]).flatten()
predictions[fn] = iffo.dot(y_pred_common_i)
return predictions