Contact
Un de nos experts vous contactera dans les plus brefs délais pour discuter de vos besoins et organiser une démo sur-mesure de la plateforme.
Les modèles d'ordre réduit (ROM) sont des modèles légers en calcul dérivés de modèles physiques plus complexes. Le but de cette réduction est de diminuer le temps nécessaire à la simulation des modèles. Ceci est généralement réalisé en résolvant un problème donné plusieurs fois sous un ensemble de paramètres, puis en déduisant de cela un modèle plus simple.
Les ROM peuvent avoir différentes utilisations :
Comme vous pouvez déjà l'imaginer, étant donné que cette technique implique d'effectuer un certain nombre de simulations d'entraînement, le calcul distribué dans le cloud est très pertinent dans ce cas.
Dans cet article, nous verrons comment fonctionne la principale technique de réduction. Ensuite, nous présenterons une application concrète de cette technique et une implémentation basée sur Python et les bibliothèques PyMOR et FeniCS, avec des calculs effectués sur la plateforme de Qarnot.
Dans cette section, nous présenterons la principale technique de réduction. Cela nécessitera d'expliquer comment les équations physiques sont généralement résolues par des techniques numériques. Ne vous inquiétez pas si vous n'êtes pas familier avec les mathématiques, nous ne plongerons pas dans les détails ici, mais présenterons uniquement les idées utilisées.
Définissons d'abord le type de problème que nous essayons de résoudre. Les phénomènes physiques tels que la propagation de la chaleur à travers une pièce métallique, l'écoulement de l'air autour d'un avion ou la propagation des ondes électromagnétiques sont tous décrits mathématiquement par des équations. Celles-ci décrivent la relation entre les variations dans le temps et l'espace des quantités physiques. Par exemple, l'équation de la chaleur nous dit comment la température change dans le temps en fonction de sa distribution spatiale à un instant donné. En trois dimensions, elle peut s'écrire :
ρcp∂t∂T=∇⋅(λ∇T)+p
Il y a deux choses à garder à l'esprit à propos de ces équations :
Il est généralement impossible de trouver mathématiquement la solution d'une telle équation aux dérivées partielles. En fait, certaines EDP courantes sont si difficiles que les mathématiciens ont du mal à démontrer quoi que ce soit à leur sujet. C'est le cas de l'équation de Navier-Stokes (décrivant l'écoulement d'un fluide) où il n'a pas encore été démontré qu'une solution régulière existe toujours (un prix d'1 million est d'ailleurs en jeu sur ce problème).
Par conséquent, les ingénieurs et les scientifiques recherchent des solutions approchées par des moyens numériques. Dans ces problèmes, certains paramètres peuvent entrer en jeu et peuvent changer. Il peut s'agir par exemple de propriétés physiques du matériau qui peuvent changer si nous décidons de remplacer un matériau par un autre ou si nous étudions l'utilisation de différentes possibilités. Les paramètres peuvent également être la température ambiante, ou la vitesse de l'avion, etc. Pour ces raisons, nous pouvons avoir la même géométrie et le même modèle, mais nous aurions besoin de calculer les solutions chaque fois que nous voulons étudier ce qui se passe pour un nouvel ensemble de paramètres.
Pour résoudre les équations numériquement, la première étape est toujours de dériver un modèle fini à partir de l'équation infinie donnée. En effet, comme nous l'avons vu, les modèles que nous essayons de résoudre contiennent un nombre infini d'inconnues puisque nous devons trouver les valeurs à chaque point et à chaque instant. Cela ne peut pas être fait par un ordinateur, car il ne dispose que d'une mémoire et d'une puissance de calcul finies. Par conséquent, nous devons rendre ce problème fini. Plusieurs techniques peuvent être utilisées pour effectuer ce processus de discrétisation, mais elles partagent souvent les deux mêmes étapes : le maillage et la discrétisation de l'équation. Ces deux étapes sont expliquées ci-dessous.
L'idée de cette étape est de dire que puisque nous ne pouvons pas calculer les valeurs à chaque point et à chaque instant, nous ne les calculerons qu'à certains points et instants. Ensuite, lorsque nous aurons besoin de connaître la valeur à un point donné, nous l'inférerons en utilisant les points voisins où les valeurs ont été calculées. Par exemple, elle peut être calculée par interpolation linéaire (par exemple, en prenant la moyenne pondérée des points voisins).
Un exemple de base est le suivant. Considérez une mince barre de métal chauffée de longueur 10. La pièce n'a alors qu'une seule dimension, et nous pourrions diviser la longueur en un certain nombre de segments (disons 5). Si nous connaissons la température à chacun des points d'extrémité, nous pouvons créer une courbe linéaire qui sera proche de la solution réelle.
L'idée centrale de ce processus est que maintenant, les solutions que nous recherchons sont entièrement décrites par un nombre fini d'informations. Ici, ce sont les valeurs à 0, 2, 4, 6, 8 et 10. Ainsi, nous pouvons dire que l'espace dans lequel nous recherchons une solution a maintenant une dimension finie (6 ici), ou que la solution a un nombre fini de degrés de liberté (DDL).
Dans les dimensions supérieures, la discrétisation de l'espace implique souvent de construire un maillage, qui est une partition de l'objet en petits triangles ou tétraèdres. Sur chacune de ces primitives, la fonction finale sera linéaire.
D'autre part, le temps est très facilement discrétisé en coupant la fenêtre de temps totale en petits pas de temps.
Pour comprendre les techniques d'ordre réduit, il est important de noter qu'il est analogue de dire que la fonction est une interpolation des valeurs à certains points et qu'elle est une combinaison linéaire de quelques fonctions simples. Si nous reprenons notre exemple, nous pouvons voir que la solution totale est la somme des six fonctions représentées en points :
Chacune de ces fonctions de base est la version mise à l'échelle d'une fonction plus « primitive » dont les valeurs sont 1 à un point et 0 à tous les autres points (ici la valeur est 10, pour que nous puissions mieux voir).
Dans le langage mathématique, ces 6 fonctions de base sont appelées une base de notre espace de solutions. Cela signifie que nous recherchons des solutions qui peuvent être exprimées comme des combinaisons linéaires (c'est-à-dire une somme de ces fonctions multipliées chacune par un facteur d'échelle) de ces fonctions.
Maintenant que notre problème a été restreint à un nombre limité de degrés de liberté, nous devons adapter l'équation à ce type de solutions. Il existe un principe en mathématiques qui stipule que pour qu'un système ait une solution unique, il doit y avoir autant d'équations que de degrés de liberté.
Différentes techniques peuvent être utilisées pour ce faire. La plus simple consisterait à réécrire l'équation à chaque point, en approximant la dérivée spatiale par (T(x+Δx)−T(x))/Δx, c'est-à-dire par la différence avec les points voisins. C'est ce qu'on appelle la méthode des différences finies et bien qu'elle soit la plus simple, elle n'est que rarement utilisée. Les plus utilisées sont la méthode des éléments finis (MEF) et la méthode des volumes finis (principalement pour la dynamique des fluides).
Vous trouverez ci-dessous une brève explication de la méthode des éléments finis. Si les détails ne vous intéressent pas, ou si vous la connaissez déjà, vous pouvez passer à la section suivante qui explique l'idée centrale de la réduction de l'ordre du modèle.
Pour comprendre l'idée derrière la méthode des éléments finis, nous pouvons revenir à notre équation de la chaleur pour une barre métallique et à la technique plus naturelle des différences finies. Dans cette technique, nous approximons la dérivée en un point par la différence avec les points voisins : ∂T/∂x=(T(x+Δx)−T(x))/Δx où Δx est la distance avec le point suivant. Le problème est que, comme nous le voyons, nous divisons notre différence par Δx pour avoir une pente. Si l'on veut un maillage fin pour un calcul précis, alors Δx sera très petit. Parce que nous divisons par ce très petit nombre, ∂T/∂x peut devenir très grand et très sensible aux erreurs numériques. Dans l'équation de la chaleur, c'est encore pire car nous considérons la dérivée seconde, qui est en (Δx)2.
Par conséquent, la méthode des différences finies est assez instable. Un de ses autres inconvénients est qu'il n'y a pas beaucoup de théorèmes sur cette technique. Par exemple, nous avons peu de résultats sur l'erreur que nous faisons en utilisant cette technique.
Pour lutter contre le problème d'instabilité des différences finies, l'idée de la méthode des éléments finis est d'utiliser l'opération opposée à la dérivation, l'intégration. Puisque l'intégration opère comme une moyenne, il semble logique qu'elle conduise à des modèles plus stables. Cependant, l'intégration sur l'ensemble des domaines entraîne une perte d'information. En effet, ce n'est pas parce que deux fonctions ont la même moyenne qu'elles sont égales. Pour résoudre ce problème et préserver la richesse de notre équation de fonction, nous multiplions les deux côtés par une fonction, appelée fonction test (écrite v dans la plupart des cas). Elle peut être vue comme une fonction de fenêtrage qui ne conservera qu'une partie de la fonction inconnue. Cette fois, si deux équations ont la même moyenne sur chaque fenêtre, on peut s'attendre à ce qu'elles soient égales. En fait, sous les bonnes hypothèses, nous avons une équivalence entre les deux formulations. Par exemple, en considérant l'équation de la chaleur en régime permanent, les deux formulations suivantes sont équivalentes (où V est l'espace approprié des fonctions test) :
∇⋅(λ∇T)+p=0 est équivalent àˋ∫Ω(∇⋅(λ∇T)+p)vdV=0,∀v∈V
L'équation de droite est appelée l'équation sous forme faible (weak form equation) et est celle qui est utilisée dans la méthode des éléments finis. Après avoir modifié l'expression, l'idée clé de la MEF n'est pas seulement de chercher une solution dans un espace Vh de dimension finie n, mais aussi de la tester avec une fonction v dans ce même espace. Étant donné que Vh est de dimension finie et que l'équation est linéaire en v, vérifier l'équation pour tout v dans Vh équivaut à la vérifier pour tous les éléments de sa base (v1,…,vn). Cela conduit naturellement à un ensemble de n équations. En considérant les n coefficients λ1,…,λn de sorte que T=∑λivi, nous obtenons un système de n équations avec n inconnues. Ce système peut être résolu de manière systématique.
En conclusion, la méthode des éléments finis consiste à réécrire l'équation sous une autre forme, la forme faible, qui implique une fonction test. La discrétisation est effectuée en utilisant des fonctions test qui ont les mêmes degrés de liberté que l'espace de solution. Pour des raisons mathématiques, il se trouve que, parce que nous avons utilisé une formulation faible, nous avons beaucoup plus de résultats mathématiques et une bien meilleure compréhension de cette technique. Par exemple, le lemme de Céa donne une limite supérieure pour l'erreur que nous faisons en utilisant cette technique de discrétisation. Plus précisément, il donne une borne supérieure pour la distance entre la solution réelle et la meilleure approximation de cette solution dans notre espace de solutions.
L'idée centrale de la réduction de l'ordre du modèle est de rechercher une solution dans une base spécifique qui ne contient que quelques éléments, mais dont nous savons qu'elle suffira à décrire les solutions. Le type de base illustré dans la section sur le maillage est très pertinent lorsque nous ne savons pas à quoi ressemblera la solution. En effet, il offre une bonne flexibilité et nous sommes sûrs que nous serons en mesure de bien décrire n'importe quelle fonction sur le domaine. Pourtant, il est aussi assez naïf :
.jpeg.Ainsi, l'idée est de simuler certaines solutions en utilisant un modèle d'ordre complet (Full Order Model - FOM), puis d'en déduire une base plus petite qui sera capable d'exprimer des solutions proches de celles d'ordre complet. Ensuite, nous n'aurons qu'à dériver des équations dans ce nouvel espace de solutions. Étant donné que cet espace a un nombre limité de degrés de liberté, le modèle dérivé ne consistera qu'en ce même faible nombre d'équations. Par conséquent, il sera presque immédiat de le résoudre. De plus, puisque l'espace de solution a été construit pour pouvoir décrire des solutions proches de la solution réelle, nous savons que la solution dérivée de notre modèle d'ordre réduit sera proche de celle d'ordre complet.
Il n'y a pas beaucoup plus à savoir sur les bases de la réduction de l'ordre du modèle. Nous pouvons maintenant explorer quelques façons de construire la base de réduction. La façon la plus simple de la construire est simplement de faire quelques simulations et d'utiliser directement les solutions pour former la base (il suffirait d'orthonormaliser la base). Cependant, il est possible de faire mieux. Vous trouverez ci-dessous une brève description de certains algorithmes qui peuvent être utilisés.
Cette technique est proche d'une technique appelée ACP (analyse en composantes principales), qui est une technique très courante en science des données. L'idée est d'extraire une base de dimension inférieure qui décrit le mieux possible certains vecteurs. Ainsi, lors de l'utilisation du POD, le pipeline consiste d'abord à effectuer un certain nombre de simulations (disons 100), puis à extraire une base plus petite (disons 20) qui décrit le mieux possible les résultats obtenus.
Dans cet algorithme, l'idée est de construire la base étape par étape. En partant d'un ensemble de résultats, nous en choisirons un qui sera notre premier vecteur de la base de réduction. Le modèle d'ordre réduit lié à la base est créé et l'on peut calculer les solutions pour toutes les autres valeurs de l'ensemble de paramètres. Ensuite, la solution avec la plus grande différence entre l'ordre complet et l'ordre réduit est ajoutée à la base. Ce faisant, nous diminuons l'erreur pour la solution qui était la moins bien décrite. Cette opération est ensuite répétée jusqu'à ce que nous atteignions le nombre de vecteurs souhaité.
Les deux algorithmes susmentionnés nécessitaient de calculer les solutions pour chaque paramètre d'un ensemble d'entraînement. Dans l'algorithme glouton faible, nous suivrons l'approche de l'algorithme glouton fort, mais nous n'aurons pas besoin de calculer toutes les solutions exactes car l'erreur n'est pas calculée, mais seulement estimée à partir du ROM.
Pour en savoir plus sur ces algorithmes, je vous recommande de lire ce tutoriel, qui explique plus en détail leur fonctionnement, leurs avantages et leurs inconvénients.
Dans cette section, nous présentons maintenant une application concrète de cette technique où nous construisons efficacement un modèle réduit pour une simulation thermique, qui pourrait être utilisé ultérieurement pour l'optimisation de la conception. Par exemple, grâce à ce modèle, nous pourrons émuler très rapidement un changement de matériau.
Le problème que nous allons étudier est le refroidissement d'une puce électronique. Le domaine est composé d'une puce chauffante logée dans un boîtier qui touche un dissipateur thermique, séparé uniquement par une fine couche de pâte thermique. Nous ne simulerons pas le flux d'air, mais considérerons uniquement une loi de refroidissement de Newton. Elle stipule que le flux de chaleur entre la surface d'un solide et un fluide pour une petite section dS est dq=−h(T−Text)dS. Ici, Text est la température de l'air.
La géométrie est la suivante :
Les paramètres considérés sont :
Le modèle est résolu dans le temps pour une durée de 10 minutes, ce qui est suffisant pour atteindre l'état stationnaire.
L'ensemble du pipeline réalisé dans cette implémentation utilise 3 outils :
La première étape de la construction du modèle d'ordre réduit est de construire le modèle de grande dimension pour avoir des solutions échantillon. Pour ce faire, nous utilisons FeniCS. Tout le code relatif au modèle est écrit dans le fichier model.py.
La première étape de la création du modèle d'ordre complet est la création du maillage. Dans notre cas, nous le chargeons simplement à partir d'un fichier. Le maillage a été réalisé à l'aide du modeleur CAO FreeCAD en complément de ses modules FEM intégrés qui appellent un mailleur dédié. Dans le code suivant, nous commençons par importer dolfin qui est l'interface Python de FeniCS. Après avoir défini certaines constantes géométriques, nous définissons une fonction pour charger le maillage.
import dolfin as df
from pymor.models.basic import InstationaryModel
''''These are the geometric constants that were used for the mesh creation
They are needed for subdomains definition'''
die_l = 20/1000
die_h = 3/1000
die_w = die_l
casing_l = 30/1000
casing_h = 5/1000
casing_w = casing_l
paste_h = 0.2/1000
block_l = 40/1000
block_h = 5/1000
block_w = 40/1000
domain_dict = {'block':1,
'paste':2,
'casing':3,
'die':4}
def load_mesh(filename):
return df.Mesh(filename)Une fois le maillage chargé, nous devons définir les sous-domaines pour les différents matériaux et conditions aux limites. Cela se fait dans make_measures(mesh) qui renvoie deux objets appelés measures (mesures) qui contiennent des informations sur les sous-domaines.
def make_measures(mesh):
'''Defines subdomains for the mesh and returns associated measures'''
tol = 1e-7 # Constant for numerical precision
boundary_air = df.CompiledSubDomain('on_boundary && x[2] > 0 - tol',
tol=tol)
sub_die = df.CompiledSubDomain('x[0] <= xb + tol && x[0] >= -xb - tol
&& x[1] <= yb + tol && x[1] >= -yb - tol
&& x[2] <= z1 + tol && x[2] >= z0 - tol',
tol=tol, xb = die_l/2, yb = die_w/2,
z1 = -paste_h-casing_h+die_h, z0 = -paste_h-casing_h)
sub_casing = df.CompiledSubDomain('x[0] <= xb + tol && x[0] >= -xb - tol
&& x[1] <= yb + tol && x[1] >= -yb - tol
&& x[2] <= z1 + tol && x[2] >= z0 - tol',
tol=tol, xb = casing_l/2, yb = casing_w/2,
z1 = -paste_h, z0 = -paste_h-casing_h)
sub_paste = df.CompiledSubDomain('x[0] <= xb + tol && x[0] >= -xb - tol
&& x[1] <= yb + tol && x[1] >= -yb - tol
&& x[2] <= z1 + tol && x[2] >= z0 - tol',
tol=tol, xb = casing_l/2, yb = casing_w/2,
z1 = 0, z0 = -paste_h)
markers = df.MeshFunction('size_t', mesh, 3)
markers.set_all(domain_dict['block'])
sub_casing.mark(markers, domain_dict['casing'])
sub_die.mark(markers, domain_dict['die'])
sub_paste.mark(markers, domain_dict['paste'])
boundary_markers = df.MeshFunction('size_t', mesh, 2)
boundary_markers.set_all(0)
boundary_air.mark(boundary_markers, 1)
dx = df.Measure('dx', domain=mesh, subdomain_data=markers)
ds = df.Measure('ds', domain=mesh, subdomain_data=boundary_markers)
return dx, dsEnsuite, nous devons écrire l'équation de la chaleur dans FeniCS. Nous allons d'abord établir sa forme forte puis sa forme faible. Cette dernière est nécessaire pour FeniCS.
L'équation de la chaleur provient de la combinaison de deux équations. La première stipule que l'augmentation de la température est proportionnelle à la chaleur reçue qui peut être écrite comme la chaleur produite au point (écrite p ici) et la différence entre les flux de chaleur entrants et sortants (qui peut être écrite −∇⋅j où j est le flux de chaleur) :
ρcp∂t∂T=p−∇⋅j
La deuxième équation est la loi de Fourier, qui stipule que le flux de chaleur est proportionnel au gradient de température (qui peut être écrit ∇T) :
j=−λ∇T
En combinant ces deux équations, nous obtenons l'équation de la chaleur (ΔT est l'opérateur de Laplace appliqué à T qui est un raccourci pour ∇⋅(∇T)) :
ρcp∂t∂T=∇⋅(λ∇T)+p
Les conditions aux limites sont les suivantes, où Tair est la température de l'air. Γair désigne la limite du domaine où le dissipateur thermique est en contact avec l'air. L'expression vient naturellement de la loi de refroidissement de Newton.
−λ∇T⋅n=h(T−Tair)surΓair
Puisque T−Tair ne dépend pas de la température de l'air, nous pouvons simplifier avec Tair=0. Nous considérons également la condition initiale T=Tair=0. Enfin, p est défini constant égal à Φ à l'intérieur de la puce et 0 ailleurs.
Pour obtenir la forme faible, nous multiplions les deux côtés par la fonction test v et intégrons :
∫Ωρcp∂t∂TvdV=∫Ω∇⋅(λ∇T)vdV+∫ΩpvdV
Ensuite, nous utilisons l'identité de Green :
∫Ω∇⋅(λ∇T)vdV=∫Γλ∇T⋅nvdS−∫Ωλ∇T⋅∇vdV
En utilisant l'expression des conditions aux limites et la simplification Tair=0, nous obtenons la forme finale :
∫Ωρcp∂t∂TvdV+∫Ωλ∇T⋅∇vdV+∫ΓairhTvdS=∫ΩpvdV
Le premier terme contient l'inertie du système et conduit à la matrice de masse. Le second représente les phénomènes de diffusion et le troisième la condition aux limites avec l'air. Le quatrième terme est le terme source. Ici, cp, ρ, λ ne dépendent que du matériau, ils sont donc constants sur chaque sous-domaine. Cela signifie que nous pouvons écrire chaque intégrale comme une somme d'intégrales sur chaque sous-domaine. Par exemple, du côté de la masse, nous pouvons écrire :
∫Ωρcp∂t∂TvdV=i∈domain_dict∑ρicp,i∫Ωi∂t∂TvdV
Dans la fonction make_pymor_bindings, nous définissons chaque partie intégrale sur chaque sous-domaine et déterminons leur matrice grâce au code suivant :
def make_pymor_bindings(V, dx, ds, solver_options):
import pymor.basic as pmb
from pymor.bindings.fenics import FenicsVectorSpace, FenicsMatrixOperator
from pymor.algorithms.timestepping import ImplicitEulerTimeStepper
space = FenicsVectorSpace(V)
# Operators matrix
u = df.TrialFunction(V)
v = df.TestFunction(V)
mass_mat = [df.assemble(df.inner(u, v)*dx(i)) for i in domain_dict.values()]
diff_mat = [df.assemble(df.inner(df.grad(u), df.grad(v))*dx(i)) for i in domain_dict.values()]
robin_left_mat = df.assemble(u*v*ds(1))
source_mat = df.assemble(v*dx(domain_dict['die']))
# Function continues here with pymor bindings (see below)
Remarquez l'utilisation de dx(i) pour désigner une intégration sur le sous-domaine i. Maintenant, la matrice de masse totale serait quelque chose comme : mass_mat_tot = sum([c_p[i]*rho[i]*mass_mat[i] for i in domain_dict.values()]). Ensuite, nous pourrions construire la matrice de diffusion totale et résoudre le système total. Cependant, nous n'allons pas le faire directement et utiliserons plutôt les fonctionnalités de pyMOR pour gérer la partie paramétrique.
C'est pourquoi nous quittons FeniCS à ce stade (en fait, nous créons également quelques matrices pour calculer les normes d'erreur, ce qui sera utile plus tard). Les seules autres parties où FeniCS est utilisé sont à des fins de lecture et d'écriture.
Maintenant que nous avons créé les opérateurs de base, nous passons à l'interface pyMOR de niveau supérieur. Le fonctionnement de pyMOR passe par des interactions avec des objets abstraits appelés Opérateur (Operator). Ces objets encapsulent les matrices et les vecteurs utilisés par les solveurs.
Par conséquent, la première chose est de créer un objet opérateur pyMOR pour chaque opérateur FeniCS que nous avons créé précédemment (pmb est le module pymor) :
# Operators
mass_op = [FenicsMatrixOperator(M, V, V, solver_options=solver_options) for M in mass_mat]
diff_op = [FenicsMatrixOperator(M, V, V, solver_options=solver_options) for M in diff_mat]
robin_left_op = FenicsMatrixOperator(robin_left_mat, V, V, solver_options=solver_options)
source_op = pmb.VectorOperator(space.make_array([source_mat]))Maintenant, nous créons deux listes de paramètres, lamb et capa, dont les composants sont les différents λi et ρicp,i.
project_lamb = [pmb.ProjectionParameterFunctional(f'lamb_{dom}', 1, 0) for dom in domain_dict.keys()]
project_capa = [pmb.ProjectionParameterFunctional(f'capa_{dom}', 1, 0) for dom in domain_dict.keys()]Ensuite, nous pouvons créer les opérateurs finaux : mass, op (diffusion + condition aux limites de l'air), rhs (source) :
mass = pmb.LincombOperator(mass_op, project_capa)
op = pmb.LincombOperator([*diff_op,
robin_left_op],
[*project_lamb,
pmb.ExpressionParameterFunctional('h[0]', {'h':1})])
rhs = pmb.LincombOperator([source_op],
[pmb.ExpressionParameterFunctional('phi[0]', {'phi':1})])Ici, nous avons introduit deux nouveaux paramètres dans pyMOR, h et ϕ, tous deux de dimension 1. L'avantage très utile que nous avons obtenu avec cette construction est que, puisque nous n'avons introduit les paramètres que du côté de pyMOR et que nous avons écrit les opérateurs comme des combinaisons linéaires d'opérateurs non paramétriques avec un coefficient qui ne dépend que des paramètres, nous sommes dans un cas de séparation des paramètres. Grâce à cela, pyMOR pourra effectuer une décomposition en ligne/hors ligne (online/offline decomposition). Cela signifie que lorsque nous projetterons le modèle pour construire notre base réduite, nous pourrons projeter toute la partie non paramétrique, de sorte que l'assemblage de notre système sera très rapide.
Continuons avec notre modélisation. L'étape suivante consiste à créer notre modèle :
end_time = 3*60
dt = 1
fom = InstationaryParametricMassModel(
T = end_time,
initial_data=space.make_array([df.project(df.Constant(0), V).vector()]),
operator=op,
rhs=rhs,
mass=mass,
time_stepper=ImplicitEulerTimeStepper(end_time//dt),
num_values= None,
products={'l2': FenicsMatrixOperator(l2_mat, V, V),
'l2_0': FenicsMatrixOperator(l2_0_mat, V, V),
'h1': FenicsMatrixOperator(h1_mat, V, V),
'h1_0_semi': FenicsMatrixOperator(h1_0_mat, V, V)},
visualizer = None
)
return fomCette étape consiste essentiellement à passer tous nos opérateurs construits précédemment à un objet Model de pyMOR. Ici, nous n'avons pas pu utiliser l'InstationaryModel standard de pyMOR car il ne prend pas en charge la matrice de masse paramétrique. Étant donné que nous avions besoin d'une masse paramétrique pour effectuer la séparation des paramètres, nous avons simplement écrit une classe enfant qui calcule la matrice de masse avant d'appeler la méthode InstationaryModel standard. Le paramètre time_stepper est un objet qui effectuera tout le travail de discrétisation temporelle pour nous. Ici, nous avons utilisé un schéma d'Euler implicite. Maintenant que nous avons notre modèle, il nous suffit de définir une valeur de paramètre μ et d'appeler fom.solve(mu) pour avoir la solution pour ces paramètres.
Nous définissons enfin make_fom qui construit l'ensemble du modèle et make_param_space qui spécifie les limites pour chaque valeur de paramètre et renvoie un objet ParameterSpace.
def make_fom(filename, option={'solver': 'cg', 'preconditioner': 'ilu'}):
'''Returns the full order model specified from the mesh filename and solver options'''
mesh = load_mesh(filename)
dx, ds = make_measures(mesh)
V = make_space(mesh)
option = {'solver': 'cg', 'preconditioner': 'ilu'}
return make_pymor_bindings(V, dx, ds, option)
def make_param_space():
''''Define ranges for parameter values and return the associated parameter space
Return:
pymor.ParameterSpace: the parametre space'''
import pymor.basic as pmb
params = pmb.Parameters({
'lamb_die': 1,
'lamb_block': 1,
'lamb_paste': 1,
'lamb_casing': 1,
'capa_die': 1,
'capa_block': 1,
'capa_paste': 1,
'capa_casing': 1,
'h': 1,
'phi': 1})
capa_range = (10**5, 10**7)
param_range = {'lamb_die': (5, 100),
'lamb_block': (100, 500),
'lamb_paste' : (0.5, 10),
'lamb_casing': (10, 200),
'capa_die': capa_range,
'capa_block': capa_range,
'capa_paste': capa_range,
'capa_casing': capa_range,
'h' : (5, 100),
'phi': (2 / die_h / die_l / die_l, 30 / die_h/ die_l / die_l)}
return pmb.ParameterSpace(params, param_range)Si nous devions effectuer la réduction de modèle localement, ce serait aussi simple que ceci :
from model import make_fom, make_param_space
import pymor.basic as pmb
fom = make_fom('mesh.xml')
param_space = make_param_space(fom)
param_set = param_space.sample_randomly(120)
T = fom.solution_space.empty()
for mu in param_set:
T.append(fom.solve(mu))
reduction_basis = pod(U_train, product=fom.h1_0_semi_product, modes=25)
reductor = pmb.InstationaryRBReductor(fom,
RB=reduction_basis,
product=fom.h1_0_semi_product)
rom = reductor.reduce()Absolument. Voici la traduction professionnelle et fidèle au texte du dernier segment de l'article, qui détaille l'étape de parallélisation sur le cloud Qarnot.
Si nous décomposons ce code, il se compose des étapes suivantes :
VectorArrays T qui est essentiellement une liste de vecteurs.Reductor pyMOR qui est utilisé pour réduire le modèle et le renvoie en tant que rom. Maintenant, si nous devions appeler rom.solve, le résultat serait donné en quelques millisecondes, tandis que fom.solve prend plusieurs minutes sur un ordinateur portable.Cependant, faire cela sur un ordinateur personnel prendrait des heures, même pour ce modèle assez simple. C'est pourquoi nous allons utiliser la plateforme de calcul dans le cloud de Qarnot pour effectuer tous les calculs de solution en parallèle.
L'idée d'utiliser Qarnot est de démarrer une poignée d'instances sur le cloud qui se partageront toutes les différentes valeurs de paramètres. Elles pourront calculer toutes les solutions en parallèle, divisant ainsi le temps nécessaire par le nombre d'instances. Ensuite, une autre instance collectera les résultats pour construire le modèle d'ordre réduit. Après cela, il sera possible de sauvegarder le ROM et de l'utiliser dans d'autres tâches pour le tester et valider son efficacité.
Avant d'entrer dans l'architecture complète des tâches, voyons ce qui est nécessaire pour envoyer une tâche de calcul sur la plateforme de Qarnot. En général, trois choses sont nécessaires :
pymor au-dessus du conteneur FeniCS officiel. Dans le code, le conteneur est spécifié par les constantes DOCKER_REPO et DOCKER_TAG.model.py. Sur Qarnot, la gestion des données se fait par des buckets S3. Il est possible de charger un fichier du disque dans un bucket ou de télécharger le contenu d'un bucket sur le disque. Un bucket peut également être spécifié pour recevoir les résultats d'une tâche.DOCKER_CMD. Dans la plupart des cas, il s'agira d'exécuter un script python envoyé avec les fichiers de ressources. Tous les scripts utilisés se trouvent dans le répertoire input.À titre d'exemple, voici un extrait du code run.py pour créer la tâche de résolution du modèle pour 120 valeurs de paramètres sur 30 nœuds (c'est la tâche train que nous verrons ensuite) :
conn = qarnot.Connection('qarnot.conf')
# [...]
DOCKER_REPO = 'qarnotlab/pymor_fenics'
DOCKER_TAG = '2020.2.0_2019.1.0'
TRAIN_PARAM_NB = 120 # Number of parameter value for training
TRAIN_INST = 30 # Number of instance that will compute in parrallel
# [...]
input_bucket = conn.create_bucket('input')
input_bucket.sync_directory('input') # adding the content of directory ./input into ressource bucket
fom_res_bucket = conn.create_bucket('fom-results')
# [...]
train_task = conn.create_task('train', 'docker-batch', TRAIN_INST, job=job)
train_task.constants['DOCKER_REPO'] = DOCKER_REPO
train_task.constants['DOCKER_TAG'] = DOCKER_TAG
# h5repack compresse the solution files to quicken data exchanges
train_task.constants['DOCKER_CMD'] = f"'python3 main.py -n {TRAIN_PARAM_NB} -o train/u &&" +
"h5repack -f GZIP=1 train/u${INSTANCE_ID}.h5 train/u${INSTANCE_ID}_c.h5'"
train_task.resources.append(input_bucket)
train_task.results = fom_res_bucket
train_task.results_whitelist = r'_c.h5' # we only want to withdraw compressed files
# [...]
train_task.submit()Le pipeline complet pour construire le modèle d'ordre réduit et le tester contient 6 tâches. Voici comment elles fonctionnent ensemble.
Le flux de travail est le suivant :
train calcule sur un certain nombre d'instances les solutions pour le modèle d'ordre complet sur un ensemble aléatoire de paramètres.rom-build construit ensuite la base de réduction et réduit le modèle avec l'algorithme POD, grâce à la solution précédemment calculée. La base et le modèle sont ensuite enregistrés sur le bucket rom.write-param écrit ensuite un ensemble de paramètres dans un fichier. Ce seront les paramètres utilisés pour valider les performances du ROM.rom-val et fom-val lisent ensuite ce fichier et résolvent le modèle (respectivement pour le modèle réduit et le modèle d'ordre complet) pour ces paramètres. Les solutions de sortie sont ensuite stockées dans deux buckets.rom-compare compare enfin les solutions et indique l'erreur commise par le ROM par rapport au modèle d'ordre complet.En plus du bucket représenté sur le graphique, toutes les tâches reçoivent également le bucket input, qui contient tous les scripts python et les données pour construire le modèle.
L'ensemble du code pour créer ces tâches et les coordonner ensemble se trouve dans run.py. C'est assez simple. Les dépendances entre les tâches sont gérées par Qarnot en créant un job et en créant toutes les tâches dans le cadre du job.
Vous pouvez télécharger tout le code ici. Pour l'exécuter, dézippez simplement le fichier, accédez au dossier contenant run.py, copiez votre fichier de configuration qarnot dans le même dossier et exécutez python3 run.py. À tout moment, vous pouvez surveiller la progression de vos tâches sur notre plateforme. Vous trouverez également de l'aide dans le fichier readme.
Dans le cas de test, j'ai utilisé 120 solutions pour construire la base de réduction. De chacune de ces solutions, 14 pas de temps sont conservés pour former la liste des solutions. Ces pas de temps sont choisis pour mettre l'accent sur le début du temps, car la phase transitoire montre plus de changements que la partie stationnaire. 50 vecteurs sont finalement conservés pour former les bases de réduction.
Pour valider notre modèle, nous résolvons le modèle pour 50 valeurs de paramètres supplémentaires et comparons les résultats entre le modèle d'ordre complet et le modèle d'ordre réduit. La pire erreur relative que nous avons est dans la tâche rom-compare :
∣∣TFOM∣∣L2(0,T;H01(Ω))∣∣TROM−TFOM∣∣L2(0,T;H01(Ω))≤2.05×10−3
Un nuage de points de l'erreur relative pour chaque simulation à chaque pas de temps se trouve dans le bucket compare :
Par conséquent, notre ROM est très précis. Nous pouvons également voir l'évolution de l'erreur par rapport au nombre de vecteurs que nous incorporons dans la base de réduction :
Nous pouvons voir la forte diminution exponentielle de l'erreur avec la taille de la base. Ceci est cohérent avec les résultats mathématiques.
L'autre chose importante à mesurer est l'amélioration de l'efficacité temporelle. Ici, nous comparons les différents temps d'exécution nécessaires pour les différentes étapes. Cela peut être vu dans la sortie du script run.py :
Nous pouvons voir que le modèle d'ordre réduit est mille fois plus rapide, ce qui prouve son efficacité et sa pertinence pour la simulation en temps réel.
De plus, le temps d'entraînement a pris 6 heures de calcul, mais grâce à notre travail de parallélisation, il n'a duré que 15 minutes. Le temps total du processus est d'environ 30 minutes.
C'est tout pour cet article !
