1// NBPROC 10
2// ff-mpirun -np 4 DDM-Schwarz-Lap-3d.edp -glut ffglut  -n 11 -k 1  -d 1 -ns -gmres 1
3/*
4  a first true parallele example fisrt freefem++
5  Ok up to 200 proc for a Poisson equation..
6  See the Doc for full explaiantion
7
8  F Hecht Dec. 2010.
9  -------------------
10 argument:
11   -glut ffglut : to see graphicaly the process
12   -n N:  set the mesh cube split NxNxN
13   -d D:  set debug flag D must be one for mpiplot
14   -k K:  to refined by K all  elemnt
15   -ns: reomove script dump
16   -gmres 0   : use iterative schwarz algo.
17          1   :  Algo GMRES on residu of schwarz algo.
18          2   :  DDM GMRES
19          3   :  DDM GMRES with coarse grid preconditionner (Good one)
20*/
21
22load "MPICG"  load "medit"  load "metis"
23include "getARGV.idp"
24include "DDM-Schwarz-macro.idp"
25//include "AddLayer3d.idp"
26include "DDM-funcs-v2.idp"
27include "cube.idp"
28include "MPIplot.idp"
29
30searchMethod=0; // more safe seach algo (warning can be very expensive in case lot of ouside point)
31// 0 by default
32assert(version >=3.11);
33real[int] ttt(10);int ittt=0;
34macro settt {ttt[ittt++]=mpiWtime();}//
35
36
37verbosity=getARGV("-vv",0);
38int vdebug=getARGV("-d",1);
39int ksplit=getARGV("-k",1);
40int nloc = getARGV("-n",20);
41string sff=getARGV("-p","");
42int gmres=getARGV("-gmres",3);
43string tsolver = getARGV("-ts","CG");
44int nC = getARGV("-N" ,max(nloc/10,5));
45int sizeoverlaps=1; // size of overlap
46bool RAS=1;
47
48if(mpirank==0 && verbosity)
49  cout << " vdebug: " << vdebug << " kspilt "<< ksplit << " nloc "<< nloc << " sff "<< sff <<"."<< endl;
50
51string sPk="P2-3gd";
52func Pk=P2;
53int Pknbcomp=1;
54
55func bool  plotMPIall(mesh3 &Th,real[int] & u,string  cm)
56{ if(vdebug) PLOTMPIALL(mesh3,Pk, Th, u,{ cmm=cm,nbiso=4,fill=1,dim=3,value=1}); return 1;}
57
58mpiComm comm(mpiCommWorld,0,0);// trick : make a no split mpiWorld
59
60
61int ipart= mpiRank(comm); // current partition number
62
63if(ipart==0)  cout << " Final N=" << ksplit*nloc << " nloc =" << nloc << " split =" << ksplit <<  endl;
64
65int[int] l111=[1,1,1,1];
66settt
67
68
69int[int,int] LL=[[1,1],[1,1],[1,1]];
70real[int,int] BB=[[0,1],[0,1],[0,1]];
71int[int] NN=[nloc,nloc,nloc];
72int[int] NNC=[nC,nC,nC];
73settt
74mesh3 Thg=Cube(NN,BB,LL);
75mesh3 ThC=Cube(NNC,BB,LL);
76
77fespace VhC(ThC,P1); // of the coarse problem..
78
79
80BuildPartitioning(sizeoverlaps,mesh3,Thg,Thi,aThij,RAS,pii,jpart,comm,vdebug)
81
82if(ksplit>1)
83{
84for(int jp=0;jp<jpart.n;++jp)
85  aThij[jp]  = trunc(aThij[jp],1,split=ksplit);
86Thi =   trunc(Thi,1,split=ksplit);
87}
88
89BuildTransferMat(ipart,mesh3,Pk,1,[0],
90                 Thi,Whi,Whij,Thij,aThij,Usend,Vrecv,jpart,vdebug)
91
92Whi ui,vi;
93
94
95
96/* the definition of the Problem .... */
97func G=1.; /* ok  */
98func F=1.; /* ok  */
99macro grad(u) [dx(u),dy(u),dz(u)] //
100// warning for Dir. BC. the last win
101varf vPb(U,V)= int3d(Thi)(grad(U)'*grad(V)) + int3d(Thi)(F*V) + on(10,U=0)+on(1,U=G) ; //');// for emacs
102varf vPbC(U,V)= int3d(ThC)(grad(U)'*grad(V))  +on(1,U=0) ; //');// for emacs
103varf vPbon10(U,V)=on(10,U=1)+on(1,U=0);
104
105varf vPBC(U,V)=on(1,U=G);
106
107
108real[int] onG10 = vPbon10(0,Whi); // on 1
109real[int] Bi=vPb(0,Whi);
110
111
112matrix Ai = vPb(Whi,Whi,solver=tsolver);
113
114DMMDeffuncAndGlobals(Lap3,comm,jpart,Whi,Vhc,1,Ai,vPbC,onG10,Pii,Usend,Vrecv,[0])
115
116Lap3CheckUpdate();
117
118Whi u=0,v;
119
120
121u[]=vPBC(0,Whi,tgv=1);
122real eps=1e-10;
123Lap3DDMSolver(Bi,u,v,gmres,eps,vdebug)
124
125
126real errg =1,umaxg;
127{
128  real umax = u[].max,umaxg;
129  real[int] aa=[umax], bb(1);
130  mpiAllReduce(aa,bb,comm,mpiMAX);
131  errg=bb[0];
132  if(ipart==0)
133    cout << " umax global  = " << bb[0] << " Wtime = " << (ttt[ittt-1]-ttt[ittt-2])  << " s " <<  " " << Lap3kiter <<  endl;
134}
135
136Lap3Saveff(sff,eps,ksplit,nloc,sizeoverlaps);