Actual source code: ex39.c
slepc-3.21.1 2024-04-26
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: */
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: {
50: FN fexp,faux,fconst,fpol;
51: PetscScalar coeffs[2];
53: PetscFunctionBeginUser;
54: PetscCall(FNCreate(PETSC_COMM_WORLD,&fexp));
55: PetscCall(FNCreate(PETSC_COMM_WORLD,&fconst));
56: PetscCall(FNCreate(PETSC_COMM_WORLD,&faux));
57: PetscCall(FNCreate(PETSC_COMM_WORLD,&fpol));
59: PetscCall(FNSetType(fexp,FNEXP));
61: PetscCall(FNSetType(fconst,FNRATIONAL));
62: coeffs[0] = -1.0;
63: PetscCall(FNRationalSetNumerator(fconst,1,coeffs));
65: PetscCall(FNSetType(faux,FNCOMBINE));
66: PetscCall(FNCombineSetChildren(faux,FN_COMBINE_ADD,fexp,fconst));
68: PetscCall(FNSetType(fpol,FNRATIONAL));
69: coeffs[0] = 1.0; coeffs[1] = 0.0;
70: PetscCall(FNRationalSetNumerator(fpol,2,coeffs));
72: PetscCall(FNSetType(fphi,FNCOMBINE));
73: PetscCall(FNCombineSetChildren(fphi,FN_COMBINE_DIVIDE,faux,fpol));
75: PetscCall(FNDestroy(&faux));
76: PetscCall(FNDestroy(&fpol));
77: PetscCall(FNDestroy(&fconst));
78: PetscCall(FNDestroy(&fexp));
79: PetscFunctionReturn(PETSC_SUCCESS);
80: }
82: int main(int argc,char **argv)
83: {
84: Mat L;
85: Vec u,w,z,yex;
86: MFN mfnexp,mfnphi;
87: FN fexp,fphi;
88: PetscBool combine=PETSC_FALSE;
89: PetscInt i,k,Istart,Iend,n=199,steps;
90: PetscReal t,tend=1.0,deltat=0.01,nrmd,nrmu,x,h;
91: const PetscReal half=0.5;
92: PetscScalar value,c,uval,*warray;
93: const PetscScalar *uarray;
95: PetscFunctionBeginUser;
96: PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
98: PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
99: PetscCall(PetscOptionsGetReal(NULL,NULL,"-tend",&tend,NULL));
100: PetscCall(PetscOptionsGetReal(NULL,NULL,"-deltat",&deltat,NULL));
101: PetscCall(PetscOptionsGetBool(NULL,NULL,"-combine",&combine,NULL));
102: h = 1.0/(n+1.0);
103: c = (n+1)*(n+1);
105: steps = (PetscInt)(tend/deltat);
106: 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: 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)":""));
109: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
110: Build the 1-D Laplacian and various vectors
111: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
112: PetscCall(MatCreate(PETSC_COMM_WORLD,&L));
113: PetscCall(MatSetSizes(L,PETSC_DECIDE,PETSC_DECIDE,n,n));
114: PetscCall(MatSetFromOptions(L));
115: PetscCall(MatGetOwnershipRange(L,&Istart,&Iend));
116: for (i=Istart;i<Iend;i++) {
117: if (i>0) PetscCall(MatSetValue(L,i,i-1,c,INSERT_VALUES));
118: if (i<n-1) PetscCall(MatSetValue(L,i,i+1,c,INSERT_VALUES));
119: PetscCall(MatSetValue(L,i,i,-2.0*c,INSERT_VALUES));
120: }
121: PetscCall(MatAssemblyBegin(L,MAT_FINAL_ASSEMBLY));
122: PetscCall(MatAssemblyEnd(L,MAT_FINAL_ASSEMBLY));
123: PetscCall(MatCreateVecs(L,NULL,&u));
124: PetscCall(VecDuplicate(u,&yex));
125: PetscCall(VecDuplicate(u,&w));
126: PetscCall(VecDuplicate(u,&z));
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: for (i=Istart;i<Iend;i++) {
134: x = (i+1)*h;
135: value = x*(1.0-x)*PetscExpReal(tend);
136: PetscCall(VecSetValue(yex,i,value,INSERT_VALUES));
137: value = PetscAbsReal(x-half)-half;
138: PetscCall(VecSetValue(u,i,value,INSERT_VALUES));
139: }
140: PetscCall(VecAssemblyBegin(yex));
141: PetscCall(VecAssemblyBegin(u));
142: PetscCall(VecAssemblyEnd(yex));
143: PetscCall(VecAssemblyEnd(u));
144: PetscCall(VecViewFromOptions(yex,NULL,"-exact_sol"));
145: PetscCall(VecViewFromOptions(u,NULL,"-initial_cond"));
147: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
148: Create two MFN solvers, for exp() and phi_1()
149: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
150: PetscCall(MFNCreate(PETSC_COMM_WORLD,&mfnexp));
151: PetscCall(MFNSetOperator(mfnexp,L));
152: PetscCall(MFNGetFN(mfnexp,&fexp));
153: PetscCall(FNSetType(fexp,FNEXP));
154: PetscCall(FNSetScale(fexp,deltat,1.0));
155: PetscCall(MFNSetErrorIfNotConverged(mfnexp,PETSC_TRUE));
156: PetscCall(MFNSetFromOptions(mfnexp));
158: PetscCall(MFNCreate(PETSC_COMM_WORLD,&mfnphi));
159: PetscCall(MFNSetOperator(mfnphi,L));
160: PetscCall(MFNGetFN(mfnphi,&fphi));
161: if (combine) PetscCall(BuildFNPhi(fphi));
162: else {
163: PetscCall(FNSetType(fphi,FNPHI));
164: PetscCall(FNPhiSetIndex(fphi,1));
165: }
166: PetscCall(FNSetScale(fphi,deltat,1.0));
167: PetscCall(MFNSetErrorIfNotConverged(mfnphi,PETSC_TRUE));
168: PetscCall(MFNSetFromOptions(mfnphi));
170: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
171: Solve the problem with the Norsett-Euler scheme
172: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
173: t = 0.0;
174: for (k=0;k<steps;k++) {
176: /* evaluate nonlinear part */
177: PetscCall(VecGetArrayRead(u,&uarray));
178: PetscCall(VecGetArray(w,&warray));
179: for (i=Istart;i<Iend;i++) {
180: x = (i+1)*h;
181: uval = uarray[i-Istart];
182: value = x*(1.0-x)*PetscExpReal(t);
183: value = value + 2.0*PetscExpReal(t) - 1.0/(1.0+value*value);
184: value = value + 1.0/(1.0+uval*uval);
185: warray[i-Istart] = deltat*value;
186: }
187: PetscCall(VecRestoreArrayRead(u,&uarray));
188: PetscCall(VecRestoreArray(w,&warray));
189: PetscCall(MFNSolve(mfnphi,w,z));
191: /* evaluate linear part */
192: PetscCall(MFNSolve(mfnexp,u,u));
193: PetscCall(VecAXPY(u,1.0,z));
194: t = t + deltat;
196: }
197: PetscCall(VecViewFromOptions(u,NULL,"-computed_sol"));
199: /*
200: Compare with exact solution and show error norm
201: */
202: PetscCall(VecCopy(u,z));
203: PetscCall(VecAXPY(z,-1.0,yex));
204: PetscCall(VecNorm(z,NORM_2,&nrmd));
205: PetscCall(VecNorm(u,NORM_2,&nrmu));
206: PetscCall(PetscPrintf(PETSC_COMM_WORLD," The relative error at t=%g is %.4f\n\n",(double)t,(double)(nrmd/nrmu)));
208: /*
209: Free work space
210: */
211: PetscCall(MFNDestroy(&mfnexp));
212: PetscCall(MFNDestroy(&mfnphi));
213: PetscCall(MatDestroy(&L));
214: PetscCall(VecDestroy(&u));
215: PetscCall(VecDestroy(&yex));
216: PetscCall(VecDestroy(&w));
217: PetscCall(VecDestroy(&z));
218: PetscCall(SlepcFinalize());
219: return 0;
220: }
222: /*TEST
224: test:
225: suffix: 1
226: args: -n 127 -tend 0.125 -mfn_tol 1e-3 -deltat 0.025
227: timeoutfactor: 2
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
237: TEST*/