FreeFEM w ICM
From Centrum Komputerów Dużej Mocy, ICM Uniwersytet Warszawski
| FreeFEM | |
|---|---|
| Produkt: | FreeFEM |
| Producent: | Université Pierre et Marie Curie |
| Licencja: | |
| Zainstalowany na: | |
| Wersja: | |
| Email: | M.Cytowski@icm.edu.pl |
| FreeFEM | |
| FreeFEM w ICM | |
| Dokumentacja | |
| Linki | |
Contents |
Instalacja
Bardzo prosimy o kontakt na adres email pomoc@icm.edu.pl, jeśli są Państwo zainteresowani prowadzeniem obliczeń w pakiecie FreeFEM w ICM.
Uruchamianie - krok po kroku
Krok 1 - problem i kod FreeFEM++
Skorzystajmy z przykładu, który oblicza równanie Allen'a-Cahn'a na kwadracie ![[0,1]\times [0,1]](/kdmwiki-dir2/images/math/3/8/c/38c7b34812cb087f3562f9ebc984e3c8.png)
. Na brzegach określamy zerowe warunki Neumanna, a jako rozwiązanie początkowe przyjmijmy nieciągła funkcję, która na trzech wybranych kołach zawartych w dziedzinie przyjmuje wartość 1, a wszędzie indziej wartość -1. Dynamika systemu w krótkim czasie doprowadzi nas do rozwiązania, które było już prezentowane na stronie głównej FreeFEM:
Na stronie [http://www.icm.edu.pl/~sheed/research/ac.html] istenieje możliwość obejrzenia filmu z tych obliczeń.
Napisany przeze mnie kod FreeFEM++ tego programu wygląda następująco:
//Allen-Cahn equation solution with FreeFEM++, written by M.Cytowski
// modify these if necessary: simulation time and output freq.
real dt = 0.00001, T = 0.002;
real eps = 0.008;
// no need to modify anything below this line
mesh Th=square(200,200);
// initial condition (initial bulks - 3 circles)
func u0=( sqrt((x-0.5)^2+(y-0.5)^2)<0.2 || sqrt( (x-0.805)^2+(y-0.5)^2)<0.1 || sqrt( (x-0.735)^2+(y-0.69)^2)<0.08 ? 1:-1);
// functions
fespace Vh(Th, P1);
Vh u=u0, uold, psi;
//Problem definition using the discretization proposed by D.Kessler, R.H.Nochetto, A.Schmidt in
//"A posteriori error control for the Allen-Cahn problem: circumventing Gronwall's inequality"
problem AC(u,psi) = int2d(Th) ( eps*u*(1/dt)*psi ) - int2d(Th) (eps*uold*(1/dt)*psi) + int2d(Th)( eps*( dx(u)*dx(psi)+dy(u)*dy(psi) ) )
+ int2d(Th)( (1/eps)*(uold*uold*uold-uold)*psi ) + int2d(Th)( (1/eps)*(3*uold*uold-1)*u*psi)
- int2d(Th)( (1/eps)*(3*uold*uold-1)*uold*psi);
//Below are the isovalues and colormap setup proposed by Lukasz Bolikowski
// isovalues
int niso = 12;
real[int] viso(niso+1);
for(int i = 0; i <= niso; i++)
viso[i] = -1.2 + 2.4 * i / niso;
viso[0] = -10000.0;
viso[niso] = 10000.0;
viso[1] = -1.0;
viso[niso-1] = 1.0;
// colormap
real[int] hsv = [ 0.0, 1.0, 1.0, 0.3, 1.0, 1.0 ];
//Adaptation error
real error=0.1;
//Adapt mesh to initial condition
Th=adaptmesh(Th,u,err=error);
// solve for T/dt timesteps with mesh adaptation after each step
for (int tt = 0; tt < int(T/dt); tt++) {
real t = tt*dt;
int numer=100000+tt;
uold = u;
AC;
plot(u, viso=viso(0:niso), hsv=hsv, fill=true, nbiso=15, cmm="[eps="+eps+"; dt="+dt+"] step="+t, ps=+numer+"ac.ps");
plot(Th,ps=+numer+"acgrid.ps");
Th=adaptmesh(Th,u,err=error);
}
Aby dowiedzieć się coś więcej o składni FreeFEM++ w celu rozszyfrowania powyższego zapisu, zapraszam na stronę dokumentacji.
Uruchomienie
Powiedzmy, że kod zapisany powyżej znajduje się w pliku ac.edp. Normalnie uruchomienie FreeFEM++ wyglądałoby następująco:
/opt/freefem++/bin/FreeFem++-nw ac.edp
Musimy jednak pamiętać, że na komputerze halo istnieje obowiązek korzystania z systemu kolejkowego PBS. Piszemy zatem następujący skrypt (plik submit.pbs):
#!/bin/tcsh #PBS -l nodes=1:ppn=1 #PBS -l file=1gb #PBS -l mem=1gb #PBS -l walltime=00:15:00 #PBS -N ac mkdir /tmp1/$user/ mkdir /tmp1/$user/ac/ echo Job $PBS_JOBNAME started echo " " at `date` echo " " on host `hostname` echo " " will run on `cat < $PBS_NODEFILE` cd /tmp1/$user/ac/ /opt/freefem++/bin/FreeFem++-nw /home/user/$user/ac.edp echo "Finished at " `date`
Następnie wstawiamy zadanie do kolejki:
qsub submit.pbs
