Department of Scientific Computing    Institute for Numerical Simulation    University of Bonn

Documentation

Download

Programming References

Bug Reports / Suggestions

FAQ

Authors

           

---

PO.cc

// adaptive solution of a Poisson-Problem with Dirichlet boundary conditions. 
// The common solve-refine- cycle is employed to successively increase accuracy.
// The threshold parameter for the adaptive refinement criterion is 1e-4. 
//
// For the discretization 4th order Interpolets and  4th order finite differences are used. 

# include "AdaptiveData.hpp"
# include "Solver.hpp"
# include "UniformData.hpp"
# include "Function.hpp"

int debugRefine ;

#define D 2 
int main() {
 Wavelets WC("Interpolet",4) ;
 
 int i,k,BC[D][2], GBC[D][2];
 // boundary condition flags:
 // 0 means Dirichlet BC at the respective face 
 // 1 means Neumann   BC at the respective face 
 //
 // try to play around a little bit with the settings
 // at least one face must have Dirichlet b.c. to make sure the linear systems have a unique solution,
 // BiCG tends to diverge otherwise ....
 BC[0][0]=0 ; // face x=0
 BC[0][1]=1 ; //      x=1 
 BC[1][0]=1 ; //      y=0
 BC[1][1]=1 ; //      y=1

 for (i=0; i<D; i++) GBC[i][0]=GBC[i][1]=-1;

 // set parameters
 AdaptivityCriterion  R(AdaptivityCriterion::L2_THRESHOLD, 1e-4) ; // refinement criterion: |u_lt| >= 1e-4
 R.MaxLevel   =LMAX-1 ;                                            // and |l|_infty < JMAX-1
 int    LL    =10     ;                                            // uniform grid for error evaluation has mesh size 2^{-10}

 printf("Refinement: MaxLevel=%d  threshold=%e\n",R.MaxLevel, R.eps) ;

 GeneralRadialFunction<D>        FU  ;  // true solution: u(x,y)=(1e-8 + (x-3/9)^2 + (y-4/9)^2)^(1/4)
 GeneralRadialFunctionDX<D>      FDX ;  // dx u(x,y), required for Neumann boundary conditions
 GeneralRadialFunctionDY<D>      FDY ;  // dy u(x,y),      "     "     "      "         "
 LaplaceGeneralRadialFunction<D> LU  ;  // Laplacian,      "     "   r.h.s.
 
 double p[4]={0.5,1e-8, 3./9 , 4./9} ;
 FU .SetParameters(p) ;
 LU .SetParameters(p) ;
 FDX.SetParameters(p) ;
 FDY.SetParameters(p) ;

 // functions for boundary values:
 // evaluate FU or FDX or FDY at the faces x[i]=k
 SubFunction FBV[D][2] ;
 double v[4]={0,D,0,0} ;
 for (i=0; i<D; i++)
   for (k=0; k<2; k++) {
     *((Function **)v) = (BC[i][k]==0) ? (Function *)&FU : ( (i==0) ? (Function *)&FDX : (Function *)&FDY ) ;
     v[2]=i ;
     v[3]=k ; 
     FBV[i][k].SetParameters((double *)v) ;
   }

 // adaptive data structures
 AdaptiveGrid<D> G(&WC) ;
 G.SetPeriodicConditions(BC) ;

 // very coarse initial grid
 G.SetFull(WC.Level0+0)  ;

 // U : solution
 // UB: function which takes the boundary  values 
 // AC: AC(l,t) >0 if the corresponding wavelet coefficient actually fulfills the refinement criterion
 AdaptiveData<D> U(&G),UB(&G),RHS(&G),AC(&G),*Tmp0=G.tmp[0], *Tmp1=G.tmp[1], *Tmp2=G.tmp[2];

 U.SetRefine()  ; // U shall be refined
 AC.SetRefine() ; // the flag field too

 U.SetBoundaryConditions (BC ) ; // some initial values
 AC.SetBoundaryConditions(GBC) ;
 U.Set(0.0)                    ;
 AC.Set(0.0)                   ;

 // AdaptiveData for boundary values
 AdaptiveData<D-1> *BV[2][2] ; 
 for (i=0; i<D; i++)
   for (k=0; k<2; k++) BV[i][k]=new AdaptiveData<D-1>(G.FaceGrid[i][k]) ;

 // operator
 HelmholtzOperator  <AdaptiveData<D>,D> Laplace ; 
 Laplace.Tmp =G.tmp[AdaptiveGrid<D>::NUM_TMP-1] ; // auxiliary AdaptiveData required for evaluation 
 Laplace.Tmp2=G.tmp[AdaptiveGrid<D>::NUM_TMP-2] ; // Apply() automatically uses 4th order finite 
 Laplace.par[0]=0 ; // 0*U +1*dxxU +1*dyyU        // differences for 4th order Interpolets!
 Laplace.par[1]=1 ;
 Laplace.par[2]=1 ;

 Solver<AdaptiveData<D>,D> S ;
 S.eps          =1e-14 ;                             // stop if residual smaller than 1e-14 
 S.maxiter      =2000  ;                             // do 200 iterations at most
 S.Op           =&Laplace ;                          // use Laplace.Apply(.,.) for Matrix-vector multiply 
 S.print        =NULL ;
 S.prstep       =1    ;
 S.StopCriterion=StoppingCriterionAbsoluteResidual ; // use this residual

 for (i=0; i<Solver<AdaptiveData<D>,D>::NUMTMP_BICGSTAB2; i++) S.Tmp[i]=G.tmp[i] ;

 double   errL1 , errL2,  errLmax, AerrL1, AerrL2, AerrLmax ;

 int LM[2]={LL,LL},  A[2]={0,0},E[2]={1<<LM[0], 1<<LM[1]} ;
 UniformData<D> M(A,E,&WC), C(A,E,&WC);
 
 int cgit=10, ll[D], innerloop=0, maxused, dof,cgitmax=5 ;

 printf("output format for L_1 / L_2 / L_infty -errors evaluated on full grids (full) and adaptive grids (adp):\n") ;
 printf("                          L_1(full)    L_1(adp)   |   L_2(full)    L_2(adp)   | L_infty(full) L_infty(adp)| max level #DOF BiCG iter\n") ;
 
 // refinement loop
 while ( (cgit > cgitmax) && (innerloop<35) ) {

   // generate right hand side
   RHS.SetFunction(&LU) ;

   // generate (interpolants) of boundary values on all faces
   for (i=0; i<D; i++) 
     for (k=0; k<2; k++) BV[i][k]->SetFunction(&FBV[i][k]) ;

   // compute function which takes the boundary values
   UB.SetBoundaryValueFunction(BV,BC,true) ;

   // put on r.h.s.
   Laplace.Apply (&UB,Tmp0) ; 
   RHS.MMinus(Tmp0)         ;

   for (i=0; i<D; i++) RHS.ApplyOp(BC[i], &RHS, i, PROJECTION) ; // project to test space (with homogeneous b.c.)

   cgit = S.BiCGSTAB2(&U,&RHS) ; // solve
          S.BiCGSTAB2(&U,&RHS) ; // re-iteration

   for (i=0; i<D; i++) Tmp0->ApplyOp(GBC[i], (i==0) ? &U : Tmp0, i, PROJECTION) ; // transform to basisfunctions without b.c.
   UB.PPlus(Tmp0) ;  // UB is complete solution including the Dirichlet-boundary values
  
   // errors and output
   //
   // Adaptive Error Calculation
   Tmp1->SetFunction(&FU) ; // true solution
   Tmp1->MMinus(&UB) ;      // error
   AerrL1  =Tmp1->L1Norm  (Tmp2) ;
   AerrL2  =Tmp1->L2Norm  (Tmp2) ;
   AerrLmax=Tmp1->LmaxNorm(Tmp2) ;
      
   // Full Grid Error Calculation
   UB.ToUniform(&M,LM) ;
   for (i=0; i<D; i++) M.ApplyOp(GBC[i], &M, i, INVERSETRANSFORM) ;

   C.SetBoundaryConditions(GBC) ;
   C.SetFunction(&FU,false,true) ;
   C.MMinus(&M) ;

   C.Abs(&C) ;
   errL1  =C.Sum() / (E[0]*E[1]) ;
   errL2  =sqrt( C.InnerProd(&C)/(E[0]*E[1]) ) ;
   errLmax=C.Max() ;

   // info on grid
   G.LargestActiveLevel(ll)     ;
   maxused=max(ll[0],ll[1])     ;
   dof=G.Size() ;

   printf("Error(%2d,%e): %e %e | %e %e | %e %e | %2d %4d %d\n", R.MaxLevel,R.eps, 
           errL1,AerrL1, errL2,AerrL2, errLmax, AerrLmax, maxused, dof, cgit) ;
 
   // refinement
   if ((cgit> cgitmax) || (innerloop==0) ) {
              
     // compute wavelet coefficients to be used for refinement control
     Tmp0->Abs(&UB) ;
     Tmp0->Add(1.0,Tmp0, 1.0,&AC) ; //entries once marked as active will again and again be marked as active ! 
     
     // refine and save the information about the active entries
     G.Refine(Tmp0 , &R, &AC) ;  
     innerloop++ ;
   } 
 } // refinement loop
  
 UB.WriteSparse("../../Data/Test/U.adp") ;
 M.WriteUDF("../../Data/Test/U") ;

}

---