libfreefem++-dev 3.47+dfsg1-2build1 / usr / include / freefem++ / idp / DDM-funcs-v2.idp

NewMacro DMMDeffuncAndGlobals(prefix,comm,jpart,Whi,Vhc,Pknbcomp,Ai,vPbC,onG10,Pii,Usend,Vrecv,U2V)
/*
    
*/

int prefix#rgmres=0;
int prefix#kiter=-1;
/*  the global name user... */
/*
  jpart,njpart : partition 
  Usend, Vrecv : buffer
  Ai , Bi 
  rMj, sMj : matrices
  
  onG10:  tgv only on DDM s Gamma and not on Gamma
*/
/*----  for coarse  solver ... */


matrix prefix#AC,prefix#Rci,prefix#Pci;/**/
{
int[int] Sigma11=U2V;   
prefix#Pci=   interpolate(Whi,VhC,U2Vc=Sigma11);
prefix#Rci =  prefix#Pci'*Pii; 
/*Rci=   interpolate(Whi,VhC,t=1,U2Vc=Sigma11(0:Pknbcomp-1));
Pci =  Pii*Rci'; */
}

/*----End of Global Def -------------*/

/*-----------------*/
/*-----------------*/


func bool prefix#Update(real[int] &ui, real[int] &vi)
{
  for(int j=0;j<jpart.n;++j)
    Usend[j][]=sMj[j]*ui; 
   SendRecvUV(comm,jpart,Usend,Vrecv)
     vi = Pii*ui;
   for(int j=0;j<jpart.n;++j)
     vi += rMj[j]*Vrecv[j][]; 
   return true;
}
/*-----------------*/
/*-----------------*/
/*-----------------*/
 
func bool  prefix#CoarseSolve(real[int]& V,real[int]& U,mpiComm& comm)
{
    if(prefix#AC.n==0 && mpiRank(comm)==0)
      prefix#AC = vPbC(VhC,VhC,solver=sparsesolver);
   /*  solvibg the coarse probleme */
   real[int] Uc(prefix#Rci.n),Bc(Uc.n); 
   Uc= prefix#Rci*U;
   mpiReduce(Uc,Bc,processor(0,comm),mpiSUM);
   if(mpiRank(comm)==0) 
      Uc = prefix#AC^-1*Bc;
    broadcast(processor(0,comm),Uc);
   V = prefix#Pci*Uc;
}

/*-----------------*/
/*-----------------*/
/*-----------------*/

func real[int] prefix#DJ(real[int]& U)
{ 
  ++prefix#kiter;
  real[int] V(U.n); 
  V =  Ai*U;
  V = onG10 ? 0.: V;  
  return V; 
}

func real[int] prefix#PDJ(real[int]& U) /* C1*/
{ 
  real[int] V(U.n); 
  
  real[int] b= onG10 ? 0. :  U; 
  V =  Ai^-1*b;	
  prefix#Update(V,U);
  return U; 
}

func real[int] prefix#PDJC(real[int]& U) /**/ 
{ /* Precon  C1= Precon , C2  precon Coarse
   Idea : F. Nataf. 
     0 ~  (I-C1A)(I-C2A) => I ~  - C1AC2A +C1A +C2A 
     New Prec P= C1+C2 - C1AC2   = C1(I- A C2) +C2
   (  C1(I- A C2) +C2 ) Uo 
     V =  - C2*Uo
  .... */
  real[int] V(U.n); 
  prefix#CoarseSolve(V,U,comm);
  V = -V; /*  -C2*Uo  */
  U  += Ai*V; /* U =  (I-A C2) Uo */
  real[int] b= onG10 ? 0. :  U; 
  U =  Ai^-1*b;	/*  ( C1( I -A C2) Uo */
  V = U -V; /**/
  prefix#Update(V,U);
  return U; 
}

/*-----------------*/
/*-----------------*/
/*-----------------*/

func real[int] prefix#PDJC2(real[int]& U) /*  bogus ???? */
{ /* Precon  C1= precon Coarse C2  precon Precon
   Idea : F. Nataf. 
      0 ~  (I C1A)(I-C2A) => I ~  - C1AC2A +C1A +C2A 
      New Prec P= C1+C2 - C1AC2   = C1(I- A C2) +C2
      (  C1(I- A C2) +C2 ) Uo 
      V =  - C2*Uo
      .... 
      V = - C2 Uo 
      W = Uo + A V 
      V + C1 W  
      */
  real[int] V(U.n);
  real[int] b= onG10 ? 0. :  U; 
  V =  Ai^-1*b;	
  b=U; 
  
  V = -V; 
  prefix#Update(V,U);
  U  += Ai*V;  

  prefix#CoarseSolve(U,b,comm);
  V = U -V;   
  prefix#Update(V,U);
  return U; 
}
/*-----------------*/
/*-----------------*/
/*-----------------*/

 func real[int] prefix#DJ0(real[int]& U)
{ 
  ++prefix#kiter;
  real[int] V(U.n); 
  real[int] b= onG10 .* U;
  b  = onG10 ? b : Bi ;  
  V = Ai^-1*b;	
  prefix#Update(V,U);
  V -= U;
   return V; 
}

/*-----------------*/
/*-----------------*/
/*-----------------*/

func bool  prefix#CheckUpdate()
{ /* verification.....*/

  Whi defPk#Pknbcomp(u,) =Times#Pknbcomp(1),defPk#Pknbcomp(v,);
  prefix#Update(u[],v[]);
  u[]-=v[];
  if(mpirank==0) cout << "CheckUpdate  " << u[].linfty << endl;
  if(u[].linfty>1e-6) plot(u,cmm="bug ????");
  assert( u[].linfty<1e-6); 
  
  return 1; }

func bool prefix#Saveff(string sff,real epss,int ksplit,int nloc,int sizeoverlaps)
{  if(sff != "")
  {
    ofstream ff(sff+".txt",append);
    cout << " ++++  " ;
    cout  << mpirank <<"/" <<  mpisize << " k=" <<  ksplit << " n= " << nloc << " " 
           << sizeoverlaps << " it=  " << prefix#kiter  ;  
    for (int i=1; i<ittt;++i)	
      cout << " " << ttt[i]-ttt[i-1] << " ";
    cout << epss << " " << Ai.nbcoef << " " << Ai.n << endl;
    
    /*
      1 mpirank
      2 mpisize 
      3 ksplit 
      4 nloc
      5 sizeoverlaps
      6 kiter
      7 mesh & part build  
      8 build the partion 
      9 build mesh, transfere , and the fine mesh .. 
      10 build the matrix,  the trans matrix, factorizatioon 
      11 GMRES 
    */
    ff   << mpirank << " " << mpisize << " " << sPk << " " ;
    ff <<  ksplit << " " << nloc << " " << sizeoverlaps << " " << prefix#kiter  ;  
    for (int i=1; i<ittt;++i)	
      ff << " " << ttt[i]-ttt[i-1] << " ";
    ff << epss << " " << Ai.nbcoef << " " << Ai.n << " " << gmres << endl;
    
  }
  return 1;
}
//
/*-----------------*/
/*-----------------*/
/*-----------------*/
/*-----------------*/
/*-----------------*/


NewMacro prefix#DDMSolver(Bi,u,v,gmres,epss,vdebug) 
{
 settt
 int ipart=mpiRank(comm);
 if(gmres==1)
  {
   prefix#rgmres=MPIAffineGMRES(prefix#DJ0,u[],veps=epss,nbiter=300,comm=comm,dimKrylov=100,verbosity=ipart ? 0: 50);
   real[int] b= onG10 .* u[];
   b  = onG10 ? b : Bi ;
   v[] = Ai^-1*b;	
   prefix#Update(v[],u[]);
  }
else if(gmres==2)
  prefix#rgmres= MPILinearGMRES(prefix#DJ,precon=prefix#PDJ,u[],Bi,veps=epss,nbiter=1000,comm=comm,
                         dimKrylov=100,verbosity=ipart ? 0: 50);
else if(gmres==3)
  prefix#rgmres= MPILinearGMRES(prefix#DJ,precon=prefix#PDJC,u[],Bi,veps=epss,nbiter=1000,comm=comm,
                          dimKrylov=100,verbosity=ipart ? 0: 50);
else if(gmres==4)
   prefix#rgmres= MPILinearGMRES(prefix#DJ,precon=prefix#PDJC2,u[],Bi,veps=epss,nbiter=1000,comm=comm,
                          dimKrylov=100,verbosity=ipart ? 0: 50);
else /*algo Shwarz for demo */
   for(int iter=0;iter <100; ++iter)
     {
       prefix#kiter=iter;
       real[int] b= onG10 .* u[];
       b  = onG10 ? b : Bi ;
       v[] = Ai^-1*b;	
       b=v[];
       prefix#Update(v[],u[]);
       if(vdebug) plotMPIall(Thi,u[],"u-"+iter);
        b -= u[];
       
       real err = b.l1;
       real umax = u[].max;
       real[int] aa=[err,umax], bb(2);
       mpiAllReduce(aa,bb,comm,mpiMAX);
       real errg = bb[0];
       real umaxg = bb[1];
       
       if(ipart==0)
	     cout << ipart << " " << iter << " err = " << errg << " u. max  " << umaxg << endl;
       if(errg< 1e-10 ) break;
     }
     prefix#Update(u[],v[]);  
 
if(vdebug) plotMPIall(Thi,v[],"u-final");
 u[]=v[];
 settt
}
EndMacro  /*prefix#DDMSolver*/

EndMacro