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.
Dans cet article, nous allons présenter FreeFEM et la Méthode des Éléments Finis. Si vous êtes un expert de FreeFEM et de la MEF, vous pouvez sauter directement à la partie Tutoriel Qarnot de cet article.
FreeFem est une implémentation open source de la Méthode des Éléments Finis Généralisée (MEFG) qui offre un moyen simple de résoudre des Équations aux Dérivées Partielles (EDP) avec des fonctions automatisées telles que la génération de maillage facile. Le projet a débuté en 1987 sous le nom de MacFem/PCFem, bien qu'il ait été réécrit en C++ et renommé FreeFem en 1992. Il est désormais maintenu par une équipe d'une douzaine d'individus, ainsi que par la communauté FreeFem.
La Méthode des Éléments Finis (MEF), utilisée par FreeFem, est appliquée à un système physique en :
La MEF est une méthode qui résout un problème physique sur un système physique donné. Les lois physiques (élasticité, électromagnétisme,...) sont modélisées par un système continu d'équations aux dérivées partielles (EDP) qui décrit la physique de l'ensemble du système. Ces EDP sont complétées par des conditions aux limites modélisant le problème aux frontières du domaine.
Ce système continu complexe est discrétisé en petites parties formant un maillage, chaque partie étant appelée un élément fini. Le choix de la taille du maillage est essentiel car il affectera la précision, les erreurs et les besoins en puissance de calcul.
À chaque élément fini est associé un système d'équations aux dérivées partielles discrétisé (représentant le système continu, adapté aux domaines discrétisés) qui peut être résolu numériquement. Ce système peut être modélisé par une matrice. Les matrices de chaque élément fini ne peuvent pas être résolues seules car elles sont interdépendantes, mais mises ensemble, elles forment une matrice représentant le système complet qui peut être résolue à l'aide d'ordinateurs.
Le système est résolu au niveau des nœuds du maillage du système, ce qui signifie que la solution n'est pas calculée en chaque point du système, mais uniquement aux nœuds, d'où le terme "éléments finis". L'interpolation, en particulier l'interpolation polynomiale, est utilisée pour résoudre le reste du système.
La MEF était initialement utilisée pour résoudre des problèmes mécaniques tels que l'analyse de la rigidité du comportement des structures d'aéronefs, mais ses applications se retrouvent dans tous les domaines de la physique numérique et de l'ingénierie.
Plusieurs étapes sont nécessaires pour exécuter une tâche sur Qarnot, mais ne vous inquiétez pas, c'est facile ! Avant de commencer, vous devez créer un compte Qarnot, nous vous offrons 15 € de calcul sur votre abonnement.
Tout d'abord, dans un dossier Qarnot_freefem_distributed_example, créez un dossier nommé resources_freefem_qarnot et enregistrez-y le script FreeFem suivant sous le nom example.edp. Cet exemple modélise un problème de "Cavité Entraînée" en 3D (Driven Cavity), un cas de Navier-Stokes. Vous pouvez bien sûr tester d'autres exemples que vous trouverez dans la documentation FreeFem très bien rédigée.
## example.edp
<details><summary>Click to see the code</summary>
<p>
```
// run with MPI: ff-mpirun -np 4 script.edp
// NBPROC 4
load "msh3"
load "MUMPS"
////////////////////
// parameters
real ttgv=1e30;
//string ssparams="nprow=1, npcol="+mpisize;
////////////////////
real nu=0.01,dt=0.3;
real alpha=1./dt,alpha2=sqrt(alpha);
int nn=7;
mesh Th2=square(nn,nn);
fespace Vh2(Th2,P2);
Vh2 ux,uz,p2;
int[int] rup=[0,2], rdown=[0,1], rmid=[1,1,2,1,3,1,4,1];
real zmin=0,zmax=1;
mesh3 Th=buildlayers(Th2,nn,
zbound=[zmin,zmax],
// region=r1,
labelmid=rmid,
reffaceup = rup,
reffacelow = rdown);
fespace VVh(Th,[P23d,P23d,P23d,P13d]);
fespace Vh(Th,P23d);
fespace Ph(Th,P13d);
macro Grad(u) [dx(u),dy(u),dz(u)]// EOM
macro div(u1,u2,u3) (dx(u1)+dy(u2)+dz(u3)) //EOM
cout << "initilisation" << endl;
real t1,t2,t3,t4;
t1=clock();
varf vStokes([u1,u2,u3,p],[v1,v2,v3,q]) =
int3d(Th,qforder=3)( Grad(u1)'*Grad(v1) + Grad(u2)'*Grad(v2) + Grad(u3)'*Grad(v3) //' for emacs
- div(u1,u2,u3)*q - div(v1,v2,v3)*p + 1e-10*q*p )
+ on(2,u1=1.,u2=0,u3=0)
+ on(1,u1=0,u2=0,u3=0)
;
matrix A=vStokes( VVh, VVh, tgv=ttgv );
t1=clock()-t1;
t4=clock();
set(A,solver=sparsesolver );
t4=clock()-t4;
t2=clock();
real[int] b= vStokes(0,VVh);
t2=clock()-t2;
VVh [u1,u2,u3,p];
VVh [X1,X2,X3,Xp];
VVh [x1,x2,x3,xp]=[x,y,z,0];
t3=clock();
u1[]= A^-1 * b;
t3=clock()-t3;
ux= u1(x,0.5,y);
uz= u3(x,0.5,y);
p2= p(x,0.5,y);
plot([ux,uz],p2,cmm=" cut y = 0.5",wait=1);
macro XX1() (x-u1*dt)//
macro XX2() (y-u2*dt)//
macro XX3() (z-u3*dt)//
varf vNS([uu1,uu2,uu3,p],[v1,v2,v3,q]) =
int3d(Th)( alpha*(uu1*v1+uu2*v2+uu3*v3) + nu*(Grad(uu1)'*Grad(v1) + Grad(uu2)'*Grad(v2) + Grad(uu3)'*Grad(v3)) //' for emacs
- div(uu1,uu2,uu3)*q - div(v1,v2,v3)*p + 1e-10*q*p )
+ on(2,uu1=1,uu2=0,uu3=0)
+ on(1,uu1=0,uu2=0,uu3=0)
+ int3d(Th,optimize=1,qforder=4)( alpha*( convect([u1,u2,u3],-dt,u1)*v1 + convect([u1,u2,u3],-dt,u2)*v2 + convect([u1,u2,u3],-dt,u3)*v3 ) ) ;
// + int3d(Th,optimize=1)( alpha*( u1(X1,X2,X3)*v1 + u2(X1,X2,X3)*v2 + u3(X1,X2,X3)*v3 ) ) ;
// + int3d(Th,optimize=1)( alpha*( u1(XX1,XX2,XX3)*v1 + u2(XX1,XX2,XX3)*v2 + u3(XX1,XX2,XX3)*v3 ) ) ;
//+ int3d(Th,optimize=1)( alpha*( u1(x,y,z)*v1 + u2(x,y,z)*v2 + u3(x,y,z)*v3 ) ) ;
//+ int3d(Th,optimize=1)( alpha*( u1*v1 + u2*v2 + u3*v3 ) ) ;
real time1,time2=0.,time3=0.,time4,time5;
time1=clock();
A = vNS( VVh, VVh, tgv=ttgv);
time1=clock()-time1;
time5=clock();
set(A,solver=sparsesolver,verb=1);
time5=clock()-time5;
time4=clock();
real t=0;
int nbtemps=5;
for(int i=0;i<nbtemps;++i)
{
real time2tmp,time3tmp;
t += dt;
cout << " iteration " << i << " t = " << t << endl;
X1[]=x1[]+u1[]*(-dt);
// verbosity=1000;
time2tmp=clock();
b=vNS(0,VVh);
time2=time2+(clock()-time2tmp);
time3tmp=clock();
u1[]= A^-1 * b;
time3=time3+(clock()-time3tmp);
ux= u1(x,0.5,y);
uz= u3(x,0.5,y);
p2= p(x,0.5,y);
plot([ux,uz],p2,cmm=" cut y = 0.5, time ="+t,wait=0);
if(i%5==6)
{
exec("mkdir dd");
string prefu="dd/pastix-nn-"+nn+"+u-"+(100+i);
string prefp="dd/pastix-nn-"+nn+"p-"+(100+i);
savemesh(Th,prefu+".mesh");
savemesh(Th,prefp+".mesh");
ofstream file(prefu+".bb");
ofstream filep(prefp+".bb");
Ph up1=u1,up2=u2,up3=u3,pp=p;
file << "3 1 3 "<< up1[].n << " 2 \n";
filep << "3 1 1 "<< pp[].n << " 2 \n";
for (int j=0;j<up1[].n ; j++)
{
file << up1[][j] <<" " <<up2[][j] <<" "<< up3[][j] <<"\n";
filep << pp[][j] << endl;
}
}
}
time4=clock()-time4;
if(mpirank == 0){
plot([ux,uz],cmm=" cut y = 0.5, time ="+t,ps="result.eps");
}
cout << "============ CPU TIME ==========================" << endl;
cout << "= Intialisation Stokes ::: =" << endl;
cout << "= matrix =" << t1<< endl;
cout << "= Factorization =" << t4 << endl;
cout << "= second member =" << t2 << endl;
cout << "= solving =" << t3<< endl;
cout << "= all =" << t1+t2+t3+t4 << endl;
cout << "============ CPU TIME ==========================" << endl;
cout << "= Navier Stokes ::: =" << endl;
cout << "= matrix =" << time1<< endl;
cout << "= Factorization =" << time5 << endl;
cout << "= second member by step (dt) =" << time2/nbtemps << endl;
cout << "= solving by step (dt) =" << time3/nbtemps<< endl;
cout << "= all step solving =" << time4 << endl;
cout << "= solving Navier Stokes =" << time1+time2+time3+time4+time5 << endl;
cout << "============ CPU TIME ==========================" << endl;
```
</p>
</details>Dans le dossier Qarnot_freefem_distributed_example, suivez ces étapes pour configurer un environnement virtuel Python. Ensuite, vous pouvez exécuter le script Python depuis votre terminal en tapant chmod +x run.py puis ./run.py. Vous pourrez alors consulter les détails des tâches sur votre propre console ou sur la console Qarnot. Une fois la tâche terminée, les résultats (un fichier de visualisation 2D) seront téléchargés sur votre ordinateur.
Nous espérons que vous avez apprécié ce tutoriel ! Si vous avez des questions ou si vous souhaitez utiliser notre plateforme pour des calculs plus lourds (nous pouvons fournir des ressources de pointe sur demande), n'hésitez pas à nous contacter.
