/usr/share/doc/freefem++/examples/examples++-tutorial/convect-apt.edp is in freefem++-doc 3.19.1-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 | // This a the rotating hill problem with one turn.
// First 1/2 turn is a convection equation and second 1/2 a convection diffusion
int kt=6;
int kloop=5;
int nbadap=5;
bool inq=0;
real tol=0.05;
real tol2=1e-4;
border a(t=0, 2*pi) { x = cos(t); y = sin(t); }; // the unit circle
mesh th = buildmesh(a(70)); // triangulates the disk
fespace Vh(th,P1);
Vh u1 = y, u2 = -x; // rotation velocity
Vh m11=0,m22=0,m12=0; // to store the metric field
Vh vT; // to save the initial in big loop
int i;
Vh vv,vo,vp; // working Finite element function
func rhill=sqrt((x-0.3)^2 +(y-0.3)^2);
func hill = 1-tanh(30*(rhill-0.2));
//func real hill(real r2) { return 1-tanh(100*(r2-0.04));}
Vh v = hill;//((x-0.3)^2 +(y-0.3)^2); // initial condition
plot(v);
vp=0;
for (int i=0;i<nbadap;i++)
{
th=adaptmesh(th,v,err=tol,inquire=inq);
v = hill;//((x-0.3)^2 +(y-0.3)^2);
real errl2= sqrt(int2d(th)(square(vp-v))) ;
vp=vp-v;
cout << " interation adp =" << i << " err l2 = " << errl2
<< " diff min = " << vp[].min << " max =" << vp[].max << endl
<< " --------------- " << endl;
vp=v;
if (errl2 < tol2) break;
}
plot(th,wait=1);
real dt = 0.17,t=0; // time step
vT[]=v[];
real error=0;
for ( i=0; i< 20/kt ; i++) {
real T=t;
vp=0;
for (int k=0;k<kloop;k++)
{
t=T; // restart
v=vT; // interpolation
m11=0;m22=0;m12=0; // reset metric
adaptmesh(th,v,err=tol,metric=[m11[],m12[],m22[]],nomeshgeneration= true ); // warning change the order in version 1.28
cout << "m11 = " << m11[].min << " " << m11[].max << endl;
cout << "m22 = " << m22[].min << " " << m22[].max << endl;
cout << "m12 = " << m12[].min << " " << m12[].max << endl;
for (int j=0;j<kt;j++)
{
t = t+dt;
vo[]=v[];
v=convect([u1,u2],-dt,vo); // convect v by u1,u2, dt seconds, results in f
plot(v,cmm="convection: t="+t + ", min=" + v[].min + ", max=" + v[].max,wait=1,dim=3);
adaptmesh(th,v,err=tol,metric=[m11[],m12[],m22[]],nomeshgeneration= true);
}
th=adaptmesh(th,v,err=tol,metric=[m11[],m12[],m22[]]);
vo=0; //
plot(th,wait=1,cmm="k=" + k + ", t= " + t + " ,i= " + i );
real errl2= sqrt(int2d(th)(square(vp-v))) ;
cout << " interation " << k << " err l2 = " << errl2 << " --------------- " <<endl;
vp=v;
error=errl2;
if (errl2 < tol2) break;
}
vT=v;
};
plot(th,v,wait=1);
|