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);