1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain

  6:    This file is part of SLEPc.
  7:    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
  8:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  9: */

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

 18: #include <slepcmfn.h>

 20: int main(int argc,char **argv)
 21: {
 22:   Mat                A;           /* problem matrix */
 23:   MFN                mfn;
 24:   FN                 f;
 25:   PetscInt           i,j,Istart,Iend,II,m,n=10,N,steps=5,its,totits=0,ncv,maxit;
 26:   PetscReal          tol,norm,h,h2,peclet=0.5,epsilon=1.0,c,i1h,j1h;
 27:   PetscScalar        t=1e-4,sone=1.0,value,upper,diag,lower;
 28:   Vec                v;
 29:   MFNConvergedReason reason;

 31:   PetscFunctionBeginUser;
 32:   PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));

 34:   PetscCall(PetscOptionsGetScalar(NULL,NULL,"-t",&t,NULL));
 35:   PetscCall(PetscOptionsGetReal(NULL,NULL,"-peclet",&peclet,NULL));
 36:   PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
 37:   PetscCall(PetscOptionsGetInt(NULL,NULL,"-steps",&steps,NULL));
 38:   m = n;
 39:   N = m*n;
 40:   /* interval [0,1], homogeneous Dirichlet boundary conditions */
 41:   h = 1.0/(n+1.0);
 42:   h2 = h*h;
 43:   c = 2.0*epsilon*peclet/h;
 44:   upper = epsilon/h2+c/(2.0*h);
 45:   diag = 2.0*(-2.0*epsilon/h2);
 46:   lower = epsilon/h2-c/(2.0*h);

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

 50:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 51:                 Generate matrix A
 52:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 53:   PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
 54:   PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N));
 55:   PetscCall(MatSetFromOptions(A));
 56:   PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
 57:   for (II=Istart;II<Iend;II++) {
 58:     i = II/n; j = II-i*n;
 59:     if (i>0) PetscCall(MatSetValue(A,II,II-n,lower,INSERT_VALUES));
 60:     if (i<m-1) PetscCall(MatSetValue(A,II,II+n,upper,INSERT_VALUES));
 61:     if (j>0) PetscCall(MatSetValue(A,II,II-1,lower,INSERT_VALUES));
 62:     if (j<n-1) PetscCall(MatSetValue(A,II,II+1,upper,INSERT_VALUES));
 63:     PetscCall(MatSetValue(A,II,II,diag,INSERT_VALUES));
 64:   }
 65:   PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
 66:   PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
 67:   PetscCall(MatCreateVecs(A,NULL,&v));

 69:   /*
 70:      Set initial condition v = 256*i^2*(1-i)^2*j^2*(1-j)^2
 71:   */
 72:   for (II=Istart;II<Iend;II++) {
 73:     i = II/n; j = II-i*n;
 74:     i1h = (i+1)*h; j1h = (j+1)*h;
 75:     value = 256.0*i1h*i1h*(1.0-i1h)*(1.0-i1h)*(j1h*j1h)*(1.0-j1h)*(1.0-j1h);
 76:     PetscCall(VecSetValue(v,i+j*n,value,INSERT_VALUES));
 77:   }
 78:   PetscCall(VecAssemblyBegin(v));
 79:   PetscCall(VecAssemblyEnd(v));

 81:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 82:                 Create the solver and set various options
 83:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 84:   PetscCall(MFNCreate(PETSC_COMM_WORLD,&mfn));
 85:   PetscCall(MFNSetOperator(mfn,A));
 86:   PetscCall(MFNGetFN(mfn,&f));
 87:   PetscCall(FNSetType(f,FNEXP));
 88:   PetscCall(FNSetScale(f,t,sone));
 89:   PetscCall(MFNSetFromOptions(mfn));

 91:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 92:                       Solve the problem, y=exp(t*A)*v
 93:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 94:   for (i=0;i<steps;i++) {
 95:     PetscCall(MFNSolve(mfn,v,v));
 96:     PetscCall(MFNGetConvergedReason(mfn,&reason));
 97:     PetscCheck(reason>=0,PETSC_COMM_WORLD,PETSC_ERR_CONV_FAILED,"Solver did not converge");
 98:     PetscCall(MFNGetIterationNumber(mfn,&its));
 99:     totits += its;
100:   }

102:   /*
103:      Optional: Get some information from the solver and display it
104:   */
105:   PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %" PetscInt_FMT "\n",totits));
106:   PetscCall(MFNGetDimensions(mfn,&ncv));
107:   PetscCall(PetscPrintf(PETSC_COMM_WORLD," Subspace dimension: %" PetscInt_FMT "\n",ncv));
108:   PetscCall(MFNGetTolerances(mfn,&tol,&maxit));
109:   PetscCall(PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4g, maxit=%" PetscInt_FMT "\n",(double)tol,maxit));
110:   PetscCall(VecNorm(v,NORM_2,&norm));
111:   PetscCall(PetscPrintf(PETSC_COMM_WORLD," Computed vector at time t=%.4g has norm %g\n",(double)PetscRealPart(t)*steps,(double)norm));

113:   /*
114:      Free work space
115:   */
116:   PetscCall(MFNDestroy(&mfn));
117:   PetscCall(MatDestroy(&A));
118:   PetscCall(VecDestroy(&v));
119:   PetscCall(SlepcFinalize());
120:   return 0;
121: }

123: /*TEST

125:    test:
126:       suffix: 1
127:       args: -mfn_tol 1e-6

129: TEST*/