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 : This example illustrates the use of Phi functions in exponential integrators.
12 : In particular, it implements the Norsett-Euler scheme of stiff order 1.
13 :
14 : The problem is the 1-D heat equation with source term
15 :
16 : y_t = y_xx + 1/(1+u^2) + psi
17 :
18 : where psi is chosen so that the exact solution is yex = x*(1-x)*exp(tend).
19 : The space domain is [0,1] and the time interval is [0,tend].
20 :
21 : [1] M. Hochbruck and A. Ostermann, "Explicit exponential Runge-Kutta
22 : methods for semilinear parabolic problems", SIAM J. Numer. Anal. 43(3),
23 : 1069-1090, 2005.
24 : */
25 :
26 : static char help[] = "Exponential integrator for the heat equation with source term.\n\n"
27 : "The command line options are:\n"
28 : " -n <idim>, where <idim> = dimension of the spatial discretization.\n"
29 : " -tend <rval>, where <rval> = real value that corresponding to the final time.\n"
30 : " -deltat <rval>, where <rval> = real value for the time increment.\n"
31 : " -combine <bool>, to represent the phi function with FNCOMBINE instead of FNPHI.\n\n";
32 :
33 : #include <slepcmfn.h>
34 :
35 : /*
36 : BuildFNPhi: builds an FNCOMBINE object representing the phi_1 function
37 :
38 : f(x) = (exp(x)-1)/x
39 :
40 : with the following tree:
41 :
42 : f(x) f(x) (combined by division)
43 : / \ p(x) = x (polynomial)
44 : a(x) p(x) a(x) (combined by addition)
45 : / \ e(x) = exp(x) (exponential)
46 : e(x) c(x) c(x) = -1 (constant)
47 : */
48 1 : PetscErrorCode BuildFNPhi(FN fphi)
49 : {
50 1 : FN fexp,faux,fconst,fpol;
51 1 : PetscScalar coeffs[2];
52 :
53 1 : PetscFunctionBeginUser;
54 1 : PetscCall(FNCreate(PETSC_COMM_WORLD,&fexp));
55 1 : PetscCall(FNCreate(PETSC_COMM_WORLD,&fconst));
56 1 : PetscCall(FNCreate(PETSC_COMM_WORLD,&faux));
57 1 : PetscCall(FNCreate(PETSC_COMM_WORLD,&fpol));
58 :
59 1 : PetscCall(FNSetType(fexp,FNEXP));
60 :
61 1 : PetscCall(FNSetType(fconst,FNRATIONAL));
62 1 : coeffs[0] = -1.0;
63 1 : PetscCall(FNRationalSetNumerator(fconst,1,coeffs));
64 :
65 1 : PetscCall(FNSetType(faux,FNCOMBINE));
66 1 : PetscCall(FNCombineSetChildren(faux,FN_COMBINE_ADD,fexp,fconst));
67 :
68 1 : PetscCall(FNSetType(fpol,FNRATIONAL));
69 1 : coeffs[0] = 1.0; coeffs[1] = 0.0;
70 1 : PetscCall(FNRationalSetNumerator(fpol,2,coeffs));
71 :
72 1 : PetscCall(FNSetType(fphi,FNCOMBINE));
73 1 : PetscCall(FNCombineSetChildren(fphi,FN_COMBINE_DIVIDE,faux,fpol));
74 :
75 1 : PetscCall(FNDestroy(&faux));
76 1 : PetscCall(FNDestroy(&fpol));
77 1 : PetscCall(FNDestroy(&fconst));
78 1 : PetscCall(FNDestroy(&fexp));
79 1 : PetscFunctionReturn(PETSC_SUCCESS);
80 : }
81 :
82 2 : int main(int argc,char **argv)
83 : {
84 2 : Mat L;
85 2 : Vec u,w,z,yex;
86 2 : MFN mfnexp,mfnphi;
87 2 : FN fexp,fphi;
88 2 : PetscBool combine=PETSC_FALSE;
89 2 : PetscInt i,k,Istart,Iend,n=199,steps;
90 2 : PetscReal t,tend=1.0,deltat=0.01,nrmd,nrmu,x,h;
91 2 : const PetscReal half=0.5;
92 2 : PetscScalar value,c,uval,*warray;
93 2 : const PetscScalar *uarray;
94 :
95 2 : PetscFunctionBeginUser;
96 2 : PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
97 :
98 2 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
99 2 : PetscCall(PetscOptionsGetReal(NULL,NULL,"-tend",&tend,NULL));
100 2 : PetscCall(PetscOptionsGetReal(NULL,NULL,"-deltat",&deltat,NULL));
101 2 : PetscCall(PetscOptionsGetBool(NULL,NULL,"-combine",&combine,NULL));
102 2 : h = 1.0/(n+1.0);
103 2 : c = (n+1)*(n+1);
104 :
105 2 : steps = (PetscInt)(tend/deltat);
106 2 : PetscCheck(PetscAbsReal(tend-steps*deltat)<10*PETSC_MACHINE_EPSILON,PETSC_COMM_WORLD,PETSC_ERR_SUP,"This example requires tend being a multiple of deltat");
107 3 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nHeat equation via phi functions, n=%" PetscInt_FMT ", tend=%g, deltat=%g%s\n\n",n,(double)tend,(double)deltat,combine?" (combine)":""));
108 :
109 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
110 : Build the 1-D Laplacian and various vectors
111 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
112 2 : PetscCall(MatCreate(PETSC_COMM_WORLD,&L));
113 2 : PetscCall(MatSetSizes(L,PETSC_DECIDE,PETSC_DECIDE,n,n));
114 2 : PetscCall(MatSetFromOptions(L));
115 2 : PetscCall(MatGetOwnershipRange(L,&Istart,&Iend));
116 256 : for (i=Istart;i<Iend;i++) {
117 254 : if (i>0) PetscCall(MatSetValue(L,i,i-1,c,INSERT_VALUES));
118 254 : if (i<n-1) PetscCall(MatSetValue(L,i,i+1,c,INSERT_VALUES));
119 254 : PetscCall(MatSetValue(L,i,i,-2.0*c,INSERT_VALUES));
120 : }
121 2 : PetscCall(MatAssemblyBegin(L,MAT_FINAL_ASSEMBLY));
122 2 : PetscCall(MatAssemblyEnd(L,MAT_FINAL_ASSEMBLY));
123 2 : PetscCall(MatCreateVecs(L,NULL,&u));
124 2 : PetscCall(VecDuplicate(u,&yex));
125 2 : PetscCall(VecDuplicate(u,&w));
126 2 : PetscCall(VecDuplicate(u,&z));
127 :
128 : /*
129 : Compute various vectors:
130 : - the exact solution yex = x*(1-x)*exp(tend)
131 : - the initial condition u = abs(x-0.5)-0.5
132 : */
133 256 : for (i=Istart;i<Iend;i++) {
134 254 : x = (i+1)*h;
135 254 : value = x*(1.0-x)*PetscExpReal(tend);
136 254 : PetscCall(VecSetValue(yex,i,value,INSERT_VALUES));
137 254 : value = PetscAbsReal(x-half)-half;
138 254 : PetscCall(VecSetValue(u,i,value,INSERT_VALUES));
139 : }
140 2 : PetscCall(VecAssemblyBegin(yex));
141 2 : PetscCall(VecAssemblyBegin(u));
142 2 : PetscCall(VecAssemblyEnd(yex));
143 2 : PetscCall(VecAssemblyEnd(u));
144 2 : PetscCall(VecViewFromOptions(yex,NULL,"-exact_sol"));
145 2 : PetscCall(VecViewFromOptions(u,NULL,"-initial_cond"));
146 :
147 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
148 : Create two MFN solvers, for exp() and phi_1()
149 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
150 2 : PetscCall(MFNCreate(PETSC_COMM_WORLD,&mfnexp));
151 2 : PetscCall(MFNSetOperator(mfnexp,L));
152 2 : PetscCall(MFNGetFN(mfnexp,&fexp));
153 2 : PetscCall(FNSetType(fexp,FNEXP));
154 2 : PetscCall(FNSetScale(fexp,deltat,1.0));
155 2 : PetscCall(MFNSetErrorIfNotConverged(mfnexp,PETSC_TRUE));
156 2 : PetscCall(MFNSetFromOptions(mfnexp));
157 :
158 2 : PetscCall(MFNCreate(PETSC_COMM_WORLD,&mfnphi));
159 2 : PetscCall(MFNSetOperator(mfnphi,L));
160 2 : PetscCall(MFNGetFN(mfnphi,&fphi));
161 2 : if (combine) PetscCall(BuildFNPhi(fphi));
162 : else {
163 1 : PetscCall(FNSetType(fphi,FNPHI));
164 1 : PetscCall(FNPhiSetIndex(fphi,1));
165 : }
166 2 : PetscCall(FNSetScale(fphi,deltat,1.0));
167 2 : PetscCall(MFNSetErrorIfNotConverged(mfnphi,PETSC_TRUE));
168 2 : PetscCall(MFNSetFromOptions(mfnphi));
169 :
170 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
171 : Solve the problem with the Norsett-Euler scheme
172 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
173 : t = 0.0;
174 12 : for (k=0;k<steps;k++) {
175 :
176 : /* evaluate nonlinear part */
177 10 : PetscCall(VecGetArrayRead(u,&uarray));
178 10 : PetscCall(VecGetArray(w,&warray));
179 1280 : for (i=Istart;i<Iend;i++) {
180 1270 : x = (i+1)*h;
181 1270 : uval = uarray[i-Istart];
182 1270 : value = x*(1.0-x)*PetscExpReal(t);
183 1270 : value = value + 2.0*PetscExpReal(t) - 1.0/(1.0+value*value);
184 1270 : value = value + 1.0/(1.0+uval*uval);
185 1270 : warray[i-Istart] = deltat*value;
186 : }
187 10 : PetscCall(VecRestoreArrayRead(u,&uarray));
188 10 : PetscCall(VecRestoreArray(w,&warray));
189 10 : PetscCall(MFNSolve(mfnphi,w,z));
190 :
191 : /* evaluate linear part */
192 10 : PetscCall(MFNSolve(mfnexp,u,u));
193 10 : PetscCall(VecAXPY(u,1.0,z));
194 10 : t = t + deltat;
195 :
196 : }
197 2 : PetscCall(VecViewFromOptions(u,NULL,"-computed_sol"));
198 :
199 : /*
200 : Compare with exact solution and show error norm
201 : */
202 2 : PetscCall(VecCopy(u,z));
203 2 : PetscCall(VecAXPY(z,-1.0,yex));
204 2 : PetscCall(VecNorm(z,NORM_2,&nrmd));
205 2 : PetscCall(VecNorm(u,NORM_2,&nrmu));
206 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," The relative error at t=%g is %.4f\n\n",(double)t,(double)(nrmd/nrmu)));
207 :
208 : /*
209 : Free work space
210 : */
211 2 : PetscCall(MFNDestroy(&mfnexp));
212 2 : PetscCall(MFNDestroy(&mfnphi));
213 2 : PetscCall(MatDestroy(&L));
214 2 : PetscCall(VecDestroy(&u));
215 2 : PetscCall(VecDestroy(&yex));
216 2 : PetscCall(VecDestroy(&w));
217 2 : PetscCall(VecDestroy(&z));
218 2 : PetscCall(SlepcFinalize());
219 : return 0;
220 : }
221 :
222 : /*TEST
223 :
224 : test:
225 : suffix: 1
226 : args: -n 127 -tend 0.125 -mfn_tol 1e-3 -deltat 0.025
227 : timeoutfactor: 2
228 :
229 : test:
230 : suffix: 2
231 : args: -n 127 -tend 0.125 -mfn_tol 1e-3 -deltat 0.025 -combine
232 : filter: sed -e "s/ (combine)//"
233 : requires: !single
234 : output_file: output/ex39_1.out
235 : timeoutfactor: 2
236 :
237 : TEST*/
|
