Line data    Source code

       1             : /*
       2             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       3             :    SLEPc - Scalable Library for Eigenvalue Problem Computations
       4             :    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
       5             : 
       6             :    This file is part of SLEPc.
       7             :    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
       8             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       9             : */
      10             : 
      11             : static char help[] = "Computes exp(t*A)*v for an advection diffusion operator with Peclet number.\n\n"
      12             :   "The command line options are:\n"
      13             :   "  -n <idim>, where <idim> = number of subdivisions of the mesh in each spatial direction.\n"
      14             :   "  -t <sval>, where <sval> = scalar value that multiplies the argument.\n"
      15             :   "  -peclet <sval>, where <sval> = Peclet value.\n"
      16             :   "  -steps <ival>, where <ival> = number of time steps.\n\n";
      17             : 
      18             : #include <slepcmfn.h>
      19             : 
      20           1 : int main(int argc,char **argv)
      21             : {
      22           1 :   Mat                A;           /* problem matrix */
      23           1 :   MFN                mfn;
      24           1 :   FN                 f;
      25           1 :   PetscInt           i,j,Istart,Iend,II,m,n=10,N,steps=5,its,totits=0,ncv,maxit;
      26           1 :   PetscReal          tol,norm,h,h2,peclet=0.5,epsilon=1.0,c,i1h,j1h;
      27           1 :   PetscScalar        t=1e-4,sone=1.0,value,upper,diag,lower;
      28           1 :   Vec                v;
      29           1 :   MFNConvergedReason reason;
      30             : 
      31           1 :   PetscFunctionBeginUser;
      32           1 :   PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
      33             : 
      34           1 :   PetscCall(PetscOptionsGetScalar(NULL,NULL,"-t",&t,NULL));
      35           1 :   PetscCall(PetscOptionsGetReal(NULL,NULL,"-peclet",&peclet,NULL));
      36           1 :   PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
      37           1 :   PetscCall(PetscOptionsGetInt(NULL,NULL,"-steps",&steps,NULL));
      38           1 :   m = n;
      39           1 :   N = m*n;
      40             :   /* interval [0,1], homogeneous Dirichlet boundary conditions */
      41           1 :   h = 1.0/(n+1.0);
      42           1 :   h2 = h*h;
      43           1 :   c = 2.0*epsilon*peclet/h;
      44           1 :   upper = epsilon/h2+c/(2.0*h);
      45           1 :   diag = 2.0*(-2.0*epsilon/h2);
      46           1 :   lower = epsilon/h2-c/(2.0*h);
      47             : 
      48           1 :   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nAdvection diffusion via y=exp(%g*A), n=%" PetscInt_FMT ", steps=%" PetscInt_FMT ", Peclet=%g\n\n",(double)PetscRealPart(t),n,steps,(double)peclet));
      49             : 
      50             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      51             :                 Generate matrix A
      52             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
      53           1 :   PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
      54           1 :   PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N));
      55           1 :   PetscCall(MatSetFromOptions(A));
      56           1 :   PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
      57         101 :   for (II=Istart;II<Iend;II++) {
      58         100 :     i = II/n; j = II-i*n;
      59         100 :     if (i>0) PetscCall(MatSetValue(A,II,II-n,lower,INSERT_VALUES));
      60         100 :     if (i<m-1) PetscCall(MatSetValue(A,II,II+n,upper,INSERT_VALUES));
      61         100 :     if (j>0) PetscCall(MatSetValue(A,II,II-1,lower,INSERT_VALUES));
      62         100 :     if (j<n-1) PetscCall(MatSetValue(A,II,II+1,upper,INSERT_VALUES));
      63         100 :     PetscCall(MatSetValue(A,II,II,diag,INSERT_VALUES));
      64             :   }
      65           1 :   PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
      66           1 :   PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
      67           1 :   PetscCall(MatCreateVecs(A,NULL,&v));
      68             : 
      69             :   /*
      70             :      Set initial condition v = 256*i^2*(1-i)^2*j^2*(1-j)^2
      71             :   */
      72         101 :   for (II=Istart;II<Iend;II++) {
      73         100 :     i = II/n; j = II-i*n;
      74         100 :     i1h = (i+1)*h; j1h = (j+1)*h;
      75         100 :     value = 256.0*i1h*i1h*(1.0-i1h)*(1.0-i1h)*(j1h*j1h)*(1.0-j1h)*(1.0-j1h);
      76         100 :     PetscCall(VecSetValue(v,i+j*n,value,INSERT_VALUES));
      77             :   }
      78           1 :   PetscCall(VecAssemblyBegin(v));
      79           1 :   PetscCall(VecAssemblyEnd(v));
      80             : 
      81             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      82             :                 Create the solver and set various options
      83             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
      84           1 :   PetscCall(MFNCreate(PETSC_COMM_WORLD,&mfn));
      85           1 :   PetscCall(MFNSetOperator(mfn,A));
      86           1 :   PetscCall(MFNGetFN(mfn,&f));
      87           1 :   PetscCall(FNSetType(f,FNEXP));
      88           1 :   PetscCall(FNSetScale(f,t,sone));
      89           1 :   PetscCall(MFNSetFromOptions(mfn));
      90             : 
      91             :   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      92             :                       Solve the problem, y=exp(t*A)*v
      93             :      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
      94           6 :   for (i=0;i<steps;i++) {
      95           5 :     PetscCall(MFNSolve(mfn,v,v));
      96           5 :     PetscCall(MFNGetConvergedReason(mfn,&reason));
      97           5 :     PetscCheck(reason>=0,PETSC_COMM_WORLD,PETSC_ERR_CONV_FAILED,"Solver did not converge");
      98           5 :     PetscCall(MFNGetIterationNumber(mfn,&its));
      99           5 :     totits += its;
     100             :   }
     101             : 
     102             :   /*
     103             :      Optional: Get some information from the solver and display it
     104             :   */
     105           1 :   PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %" PetscInt_FMT "\n",totits));
     106           1 :   PetscCall(MFNGetDimensions(mfn,&ncv));
     107           1 :   PetscCall(PetscPrintf(PETSC_COMM_WORLD," Subspace dimension: %" PetscInt_FMT "\n",ncv));
     108           1 :   PetscCall(MFNGetTolerances(mfn,&tol,&maxit));
     109           1 :   PetscCall(PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4g, maxit=%" PetscInt_FMT "\n",(double)tol,maxit));
     110           1 :   PetscCall(VecNorm(v,NORM_2,&norm));
     111           1 :   PetscCall(PetscPrintf(PETSC_COMM_WORLD," Computed vector at time t=%.4g has norm %g\n",(double)PetscRealPart(t)*steps,(double)norm));
     112             : 
     113             :   /*
     114             :      Free work space
     115             :   */
     116           1 :   PetscCall(MFNDestroy(&mfn));
     117           1 :   PetscCall(MatDestroy(&A));
     118           1 :   PetscCall(VecDestroy(&v));
     119           1 :   PetscCall(SlepcFinalize());
     120             :   return 0;
     121             : }
     122             : 
     123             : /*TEST
     124             : 
     125             :    test:
     126             :       suffix: 1
     127             :       args: -mfn_tol 1e-6
     128             : 
     129             : TEST*/