Line | Branch | Exec | Source |
---|---|---|---|
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 | 10 | int main(int argc,char **argv) | |
21 | { | ||
22 | 10 | Mat A; /* problem matrix */ | |
23 | 10 | MFN mfn; | |
24 | 10 | FN f; | |
25 | 10 | PetscInt i,j,Istart,Iend,II,m,n=10,N,steps=5,its,totits=0,ncv,maxit; | |
26 | 10 | PetscReal tol,norm,h,h2,peclet=0.5,epsilon=1.0,c,i1h,j1h; | |
27 | 10 | PetscScalar t=1e-4,sone=1.0,value,upper,diag,lower; | |
28 | 10 | Vec v; | |
29 | 10 | MFNConvergedReason reason; | |
30 | |||
31 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
10 | PetscFunctionBeginUser; |
32 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(SlepcInitialize(&argc,&argv,NULL,help)); |
33 | |||
34 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(PetscOptionsGetScalar(NULL,NULL,"-t",&t,NULL)); |
35 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(PetscOptionsGetReal(NULL,NULL,"-peclet",&peclet,NULL)); |
36 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL)); |
37 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(PetscOptionsGetInt(NULL,NULL,"-steps",&steps,NULL)); |
38 | 10 | m = n; | |
39 | 10 | N = m*n; | |
40 | /* interval [0,1], homogeneous Dirichlet boundary conditions */ | ||
41 | 10 | h = 1.0/(n+1.0); | |
42 | 10 | h2 = h*h; | |
43 | 10 | c = 2.0*epsilon*peclet/h; | |
44 | 10 | upper = epsilon/h2+c/(2.0*h); | |
45 | 10 | diag = 2.0*(-2.0*epsilon/h2); | |
46 | 10 | lower = epsilon/h2-c/(2.0*h); | |
47 | |||
48 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | 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 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(MatCreate(PETSC_COMM_WORLD,&A)); |
54 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N)); |
55 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(MatSetFromOptions(A)); |
56 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(MatGetOwnershipRange(A,&Istart,&Iend)); |
57 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
1010 | for (II=Istart;II<Iend;II++) { |
58 | 1000 | i = II/n; j = II-i*n; | |
59 |
6/8✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
|
1000 | if (i>0) PetscCall(MatSetValue(A,II,II-n,lower,INSERT_VALUES)); |
60 |
6/8✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
|
1000 | if (i<m-1) PetscCall(MatSetValue(A,II,II+n,upper,INSERT_VALUES)); |
61 |
6/8✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
|
1000 | if (j>0) PetscCall(MatSetValue(A,II,II-1,lower,INSERT_VALUES)); |
62 |
6/8✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
|
1000 | if (j<n-1) PetscCall(MatSetValue(A,II,II+1,upper,INSERT_VALUES)); |
63 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
1000 | PetscCall(MatSetValue(A,II,II,diag,INSERT_VALUES)); |
64 | } | ||
65 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); |
66 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); |
67 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | 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 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
1010 | for (II=Istart;II<Iend;II++) { |
73 | 1000 | i = II/n; j = II-i*n; | |
74 | 1000 | i1h = (i+1)*h; j1h = (j+1)*h; | |
75 | 1000 | value = 256.0*i1h*i1h*(1.0-i1h)*(1.0-i1h)*(j1h*j1h)*(1.0-j1h)*(1.0-j1h); | |
76 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
1000 | PetscCall(VecSetValue(v,i+j*n,value,INSERT_VALUES)); |
77 | } | ||
78 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(VecAssemblyBegin(v)); |
79 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(VecAssemblyEnd(v)); |
80 | |||
81 | /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - | ||
82 | Create the solver and set various options | ||
83 | - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ | ||
84 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(MFNCreate(PETSC_COMM_WORLD,&mfn)); |
85 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(MFNSetOperator(mfn,A)); |
86 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(MFNGetFN(mfn,&f)); |
87 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(FNSetType(f,FNEXP)); |
88 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(FNSetScale(f,t,sone)); |
89 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(MFNSetFromOptions(mfn)); |
90 | |||
91 | /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - | ||
92 | Solve the problem, y=exp(t*A)*v | ||
93 | - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ | ||
94 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
60 | for (i=0;i<steps;i++) { |
95 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
50 | PetscCall(MFNSolve(mfn,v,v)); |
96 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
50 | PetscCall(MFNGetConvergedReason(mfn,&reason)); |
97 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
50 | PetscCheck(reason>=0,PETSC_COMM_WORLD,PETSC_ERR_CONV_FAILED,"Solver did not converge"); |
98 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
50 | PetscCall(MFNGetIterationNumber(mfn,&its)); |
99 | 50 | totits += its; | |
100 | } | ||
101 | |||
102 | /* | ||
103 | Optional: Get some information from the solver and display it | ||
104 | */ | ||
105 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %" PetscInt_FMT "\n",totits)); |
106 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(MFNGetDimensions(mfn,&ncv)); |
107 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(PetscPrintf(PETSC_COMM_WORLD," Subspace dimension: %" PetscInt_FMT "\n",ncv)); |
108 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(MFNGetTolerances(mfn,&tol,&maxit)); |
109 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4g, maxit=%" PetscInt_FMT "\n",(double)tol,maxit)); |
110 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(VecNorm(v,NORM_2,&norm)); |
111 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | 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 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(MFNDestroy(&mfn)); |
117 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(MatDestroy(&A)); |
118 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(VecDestroy(&v)); |
119 |
3/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
|
10 | PetscCall(SlepcFinalize()); |
120 | return 0; | ||
121 | } | ||
122 | |||
123 | /*TEST | ||
124 | |||
125 | test: | ||
126 | suffix: 1 | ||
127 | args: -mfn_tol 1e-6 | ||
128 | |||
129 | TEST*/ | ||
130 |