/usr/share/doc/freefem++/examples/examples++-mpi/Stokes-v3-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 102 103 104 105 106 107 108 109 110 111 112 113 114 115 | // other
load "msh3"
load "medit"
load "real_SuperLU_DIST_FreeFem"
real ttgv=-1;
string ssparams="nprow=1, npcol="+mpisize;
int nn=12;
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;
//VVh [uu1,uu2,uu3,up];
real timeI=clock();
real time1=clock();
varf vlaplace1(u1,v1) =
int3d(Th,qforder=3)( Grad(u1)'*Grad(v1) ) + on(2,u1=fup) + on(1,u1=0); //' for emacs
varf vlaplacenull(u2,v2) =
int3d(Th,qforder=3)( Grad(u2)'*Grad(v2) ) + on(2,u2=0) + on(1,u2=0); //' for emacs
varf vdx(u1,q) = int3d(Th,qforder=3)( -dx(u1)*q );
varf vdy(u2,q) = int3d(Th,qforder=3)( -dy(u2)*q );
varf vdz(u3,q) = int3d(Th,qforder=3)( -dz(u3)*q );
varf vmass(p,q) = int3d(Th,qforder=3)( 1e-10*p*q );
matrix M = vmass(Ph,Ph,tgv=ttgv);
matrix L1 = vlaplace1(Uh,Uh,tgv=ttgv);
matrix Lnull = vlaplacenull(Uh,Uh);
matrix Bx = vdx(Uh,Ph);
matrix By = vdy(Uh,Ph);
matrix Bz = vdz(Uh,Ph);
matrix A = [ [ L1, 0, 0, Bx' ],
[ 0, Lnull, 0, By' ],
[ 0, 0, Lnull, Bz' ],
[ Bx, By, Bz, M ] ]; //' for emacs
time1=clock()-time1;
real time2=clock();
real[int] b(VVh.ndof);
real[int] xx(VVh.ndof);
{
real[int] b1 = vlaplace1(0,Uh);
real[int] b2 = vlaplacenull(0,Uh);
for(int ii=0; ii<Uh.ndof; ii++)
{
b[ii]=b1[ii];
b[ii+Uh.ndof]=b2[ii];
b[ii+2*Uh.ndof]=b2[ii];
}
for(int ii=UUh.ndof; ii<UUh.ndof+Ph.ndof; ii++)
{
b[ii]=0;
}
}
time2=clock()-time2;
real timeF=clock();
set(A,solver=sparsesolver,tgv=ttgv,sparams=ssparams);
timeF=clock()-timeF;
real time3=clock();
xx = A^-1*b;
time3=clock()-time3;
for(int ii=0; ii<Uh.ndof; ii++)
{
u1[][ii]=xx[ii];
u2[][ii]=xx[ii+Uh.ndof];
u3[][ii]=xx[ii+2*Uh.ndof];
}
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);
|