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,NULL,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*/
|
