Adaptive1DSolvers.cc

00001 //: Adaptive1D.cc
00002 
00003 #include "Adaptive1D.hpp"
00004 
00005 // Solve A*(*this)=(*RHS)
00006 int AdaptiveB::Solve(AdaptiveB *RHS, AdaptiveBOperator  *Op) {
00007   static AdaptiveB r,p,q,ap  ;
00008   double rr,rq,ak,bk,pap     ;
00009   int    it=0 ;
00010   
00011   r.Init(RHS->Level0,RHS->J) ;
00012   p.Init(RHS->Level0,RHS->J) ;
00013   q.Init(RHS->Level0,RHS->J) ;
00014   ap.Init(RHS->Level0,RHS->J) ;
00015 
00016   CopyIndexSets(RHS) ;
00017   r.CopyIndexSets(RHS) ;
00018   p.CopyIndexSets(RHS) ;
00019   q.CopyIndexSets(RHS) ;
00020   ap.CopyIndexSets(RHS) ;
00021   
00022   SetBoundaryConditions(RHS->BC) ;
00023   r.SetBoundaryConditions (RHS->BC) ;
00024   p.SetBoundaryConditions (RHS->BC) ;
00025   q.SetBoundaryConditions (RHS->BC) ;
00026   ap.SetBoundaryConditions(RHS->BC) ;
00027 
00028   Op->apply(this,&ap) ;
00029   r.Sub(RHS,&ap) ;       // r=rhs-A(*this) 
00030   
00031   Op->Preconditioner(&r,&q) ;
00032 
00033   p.Add(1.0,&q,0.0,&p) ; //p=q
00034 
00035   rr=r.InnerProd(&r) ;        // (r,r)
00036   rq=r.InnerProd(&q) ;        // (r,q)
00037 
00038   while (rr>1e-10) {
00039 
00040     Op->apply(&p,&ap) ; 
00041     pap=p.InnerProd(&ap) ;
00042 
00043     ak=rq/pap ;
00044     bk=1/rq   ;
00045 
00046     PPlus(ak,&p) ;
00047     r.PPlus(-ak,&ap)   ;
00048 
00049     Op->Preconditioner(&r,&q) ; // q=C*r 
00050 
00051     rr=r.InnerProd(&r) ;        // (r,r)
00052     rq=r.InnerProd(&q) ;    
00053 
00054     bk *= rq ;
00055     
00056     p.Add(1.0,&q , bk,&p) ;
00057     it++ ;
00058   }
00059  return it ;
00060 }

Generated at Mon Aug 19 10:02:31 2002 for AWFD by