/usr/share/doc/freefem++/examples/examples++-mpi/Stokes-v2-matrix-superludist.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 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 | load "msh3"
load "medit"
load "real_SuperLU_DIST_FreeFem"
real ttgv=1e10;
string ssparams="nprow=1, npcol="+mpisize;
int nn=10;
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], labelmid=rmid,
reffaceup = rup, reffacelow = rdown);
//medit("c10x10x10",Th);
fespace VVh(Th,[P2,P2,P2,P1]);
fespace UUh(Th,[P2,P2,P2]);
fespace Uh(Th,P2);
fespace Ph(Th,P1);
macro Grad(u) [dx(u),dy(u),dz(u)]// EOM
macro div(u1,u2,u3) (dx(u1)+dy(u2)+dz(u3)) //EOM
func fup = (1-x)*(x)*y*(1-y)*16;
Uh u1,u2,u3;
Uh v1,v2,v3;
Ph p,q;
real timeI=clock();
real time1=clock();
//VVh [uu1,uu2,uu3,up];
varf vlaplace([u1,u2,u3],[v1,v2,v3]) =
int3d(Th,qforder=3)( Grad(u1)'*Grad(v1) + Grad(u2)'*Grad(v2) + Grad(u3)'*Grad(v3) )
+ on(2,u1=fup,u2=0,u3=0) + on(1,u1=0,u2=0,u3=0) ; //'
varf vdiv([u1,u2,u3],[q])=int3d(Th,qforder=3)(- div(u1,u2,u3)*q );
varf vdivT([q],[u1,u2,u3])=int3d(Th,qforder=3)(- div(u1,u2,u3)*q );
varf vmass(p,q) = int3d(Th,qforder=3)( 1e-10*p*q );
matrix M=vmass(Ph,Ph);
matrix L=vlaplace(UUh,UUh,tgv=ttgv);
matrix B=vdiv(UUh,Ph);
matrix BT=vdivT(Ph,UUh);
matrix A = [ [ L, BT ], //'
[ B, M ]];
time1=clock()-time1;
real timeF=clock();
set(A,solver=sparsesolver,sparams=ssparams);
timeF=clock()-timeF;
real time2=clock();
real[int] b(VVh.ndof);
real[int] xx(VVh.ndof);
{
real[int] b1 = vlaplace(0,UUh);
for(int ii=0; ii<UUh.ndof; ii++)
{
b[ii]=b1[ii];
}
for(int ii=0; ii<Ph.ndof; ii++)
{
b[ii+UUh.ndof]=0;
}
}
time2=clock()-time2;
real time3=clock();
xx = A^-1*b;
time3=clock()-time3;
for(int ii=0; ii<Uh.ndof; ii++)
{
u1[][ii]=xx[3*ii];
u2[][ii]=xx[3*ii+1];
u3[][ii]=xx[3*ii+2];
}
for(int ii=0; ii<Ph.ndof; ii++)
{
p[][ii]=xx[ii+UUh.ndof];
}
timeI=clock()-timeI;
cout << "============= CPU TIME ============" << endl;
cout << " matrix " << time1 << endl;
cout << " Fact " << timeF << endl;
cout << " second member " << time2 << endl;
cout << " solve " << time3 << endl;
cout << " ------------" << endl;
cout << " all " << timeI << endl;
cout << "============= CPU TIME ============" << endl;
if(mpirank==0) medit("UV2 PV2",Th,[u1,u2,u3],p);
|