Actual source code: ex39.c
slepc-3.14.1 2020-12-08
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2020, 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: */
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.
14: The problem is the 1-D heat equation with source term
16: y_t = y_xx + 1/(1+u^2) + psi
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].
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: */
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";
33: #include <slepcmfn.h>
35: /*
36: BuildFNPhi: builds an FNCOMBINE object representing the phi_1 function
38: f(x) = (exp(x)-1)/x
40: with the following tree:
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: PetscErrorCode BuildFNPhi(FN fphi)
49: {
51: FN fexp,faux,fconst,fpol;
52: PetscScalar coeffs[2];
55: FNCreate(PETSC_COMM_WORLD,&fexp);
56: FNCreate(PETSC_COMM_WORLD,&fconst);
57: FNCreate(PETSC_COMM_WORLD,&faux);
58: FNCreate(PETSC_COMM_WORLD,&fpol);
60: FNSetType(fexp,FNEXP);
62: FNSetType(fconst,FNRATIONAL);
63: coeffs[0] = -1.0;
64: FNRationalSetNumerator(fconst,1,coeffs);
66: FNSetType(faux,FNCOMBINE);
67: FNCombineSetChildren(faux,FN_COMBINE_ADD,fexp,fconst);
69: FNSetType(fpol,FNRATIONAL);
70: coeffs[0] = 1.0; coeffs[1] = 0.0;
71: FNRationalSetNumerator(fpol,2,coeffs);
73: FNSetType(fphi,FNCOMBINE);
74: FNCombineSetChildren(fphi,FN_COMBINE_DIVIDE,faux,fpol);
76: FNDestroy(&faux);
77: FNDestroy(&fpol);
78: FNDestroy(&fconst);
79: FNDestroy(&fexp);
80: return(0);
81: }
83: int main(int argc,char **argv)
84: {
85: Mat L;
86: Vec u,w,z,yex;
87: MFN mfnexp,mfnphi;
88: FN fexp,fphi;
89: PetscBool combine=PETSC_FALSE;
90: PetscInt i,k,Istart,Iend,n=199,steps;
91: PetscReal t,tend=1.0,deltat=0.01,nrmd,nrmu,x,h;
92: const PetscReal half=0.5;
93: PetscScalar value,c,uval,*warray;
94: const PetscScalar *uarray;
95: PetscErrorCode ierr;
97: SlepcInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
99: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
100: PetscOptionsGetReal(NULL,NULL,"-tend",&tend,NULL);
101: PetscOptionsGetReal(NULL,NULL,"-deltat",&deltat,NULL);
102: PetscOptionsGetBool(NULL,NULL,"-combine",&combine,NULL);
103: h = 1.0/(n+1.0);
104: c = (n+1)*(n+1);
106: steps = (PetscInt)(tend/deltat);
107: if (PetscAbsReal(tend-steps*deltat)>10*PETSC_MACHINE_EPSILON) SETERRQ(PETSC_COMM_WORLD,1,"This example requires tend being a multiple of deltat");
108: PetscPrintf(PETSC_COMM_WORLD,"\nHeat equation via phi functions, n=%D, tend=%g, deltat=%g%s\n\n",n,(double)tend,(double)deltat,combine?" (combine)":"");
110: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
111: Build the 1-D Laplacian and various vectors
112: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
113: MatCreate(PETSC_COMM_WORLD,&L);
114: MatSetSizes(L,PETSC_DECIDE,PETSC_DECIDE,n,n);
115: MatSetFromOptions(L);
116: MatSetUp(L);
117: MatGetOwnershipRange(L,&Istart,&Iend);
118: for (i=Istart;i<Iend;i++) {
119: if (i>0) { MatSetValue(L,i,i-1,c,INSERT_VALUES); }
120: if (i<n-1) { MatSetValue(L,i,i+1,c,INSERT_VALUES); }
121: MatSetValue(L,i,i,-2.0*c,INSERT_VALUES);
122: }
123: MatAssemblyBegin(L,MAT_FINAL_ASSEMBLY);
124: MatAssemblyEnd(L,MAT_FINAL_ASSEMBLY);
125: MatCreateVecs(L,NULL,&u);
126: VecDuplicate(u,&yex);
127: VecDuplicate(u,&w);
128: VecDuplicate(u,&z);
130: /*
131: Compute various vectors:
132: - the exact solution yex = x*(1-x)*exp(tend)
133: - the initial condition u = abs(x-0.5)-0.5
134: */
135: for (i=Istart;i<Iend;i++) {
136: x = (i+1)*h;
137: value = x*(1.0-x)*PetscExpReal(tend);
138: VecSetValue(yex,i,value,INSERT_VALUES);
139: value = PetscAbsReal(x-half)-half;
140: VecSetValue(u,i,value,INSERT_VALUES);
141: }
142: VecAssemblyBegin(yex);
143: VecAssemblyBegin(u);
144: VecAssemblyEnd(yex);
145: VecAssemblyEnd(u);
146: VecViewFromOptions(yex,NULL,"-exact_sol");
147: VecViewFromOptions(u,NULL,"-initial_cond");
149: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
150: Create two MFN solvers, for exp() and phi_1()
151: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
152: MFNCreate(PETSC_COMM_WORLD,&mfnexp);
153: MFNSetOperator(mfnexp,L);
154: MFNGetFN(mfnexp,&fexp);
155: FNSetType(fexp,FNEXP);
156: FNSetScale(fexp,deltat,1.0);
157: MFNSetErrorIfNotConverged(mfnexp,PETSC_TRUE);
158: MFNSetFromOptions(mfnexp);
160: MFNCreate(PETSC_COMM_WORLD,&mfnphi);
161: MFNSetOperator(mfnphi,L);
162: MFNGetFN(mfnphi,&fphi);
163: if (combine) {
164: BuildFNPhi(fphi);
165: } else {
166: FNSetType(fphi,FNPHI);
167: FNPhiSetIndex(fphi,1);
168: }
169: FNSetScale(fphi,deltat,1.0);
170: MFNSetErrorIfNotConverged(mfnphi,PETSC_TRUE);
171: MFNSetFromOptions(mfnphi);
173: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
174: Solve the problem with the Norsett-Euler scheme
175: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
176: t = 0.0;
177: for (k=0;k<steps;k++) {
179: /* evaluate nonlinear part */
180: VecGetArrayRead(u,&uarray);
181: VecGetArray(w,&warray);
182: for (i=Istart;i<Iend;i++) {
183: x = (i+1)*h;
184: uval = uarray[i-Istart];
185: value = x*(1.0-x)*PetscExpReal(t);
186: value = value + 2.0*PetscExpReal(t) - 1.0/(1.0+value*value);
187: value = value + 1.0/(1.0+uval*uval);
188: warray[i-Istart] = deltat*value;
189: }
190: VecRestoreArrayRead(u,&uarray);
191: VecRestoreArray(w,&warray);
192: MFNSolve(mfnphi,w,z);
194: /* evaluate linear part */
195: MFNSolve(mfnexp,u,u);
196: VecAXPY(u,1.0,z);
197: t = t + deltat;
199: }
200: VecViewFromOptions(u,NULL,"-computed_sol");
202: /*
203: Compare with exact solution and show error norm
204: */
205: VecCopy(u,z);
206: VecAXPY(z,-1.0,yex);
207: VecNorm(z,NORM_2,&nrmd);
208: VecNorm(u,NORM_2,&nrmu);
209: PetscPrintf(PETSC_COMM_WORLD," The relative error at t=%g is %.4f\n\n",(double)t,(double)(nrmd/nrmu));
211: /*
212: Free work space
213: */
214: MFNDestroy(&mfnexp);
215: MFNDestroy(&mfnphi);
216: MatDestroy(&L);
217: VecDestroy(&u);
218: VecDestroy(&yex);
219: VecDestroy(&w);
220: VecDestroy(&z);
221: SlepcFinalize();
222: return ierr;
223: }
225: /*TEST
227: test:
228: suffix: 1
229: args: -n 127 -tend 0.125 -mfn_tol 1e-3 -deltat 0.025
230: timeoutfactor: 2
232: test:
233: suffix: 2
234: args: -n 127 -tend 0.125 -mfn_tol 1e-3 -deltat 0.025 -combine
235: filter: sed -e "s/ (combine)//"
236: requires: !single
237: output_file: output/ex39_1.out
238: timeoutfactor: 2
240: TEST*/