Actual source code: ex39.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2019, 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:   PetscScalar        value,c,uval,*warray;
 93:   const PetscScalar *uarray;
 94:   PetscErrorCode     ierr;

 96:   SlepcInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;

 98:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
 99:   PetscOptionsGetReal(NULL,NULL,"-tend",&tend,NULL);
100:   PetscOptionsGetReal(NULL,NULL,"-deltat",&deltat,NULL);
101:   PetscOptionsGetBool(NULL,NULL,"-combine",&combine,NULL);
102:   h = 1.0/(n+1.0);
103:   c = (n+1)*(n+1);

105:   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)":"");

107:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
108:                  Build the 1-D Laplacian and various vectors
109:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
110:   MatCreate(PETSC_COMM_WORLD,&L);
111:   MatSetSizes(L,PETSC_DECIDE,PETSC_DECIDE,n,n);
112:   MatSetFromOptions(L);
113:   MatSetUp(L);
114:   MatGetOwnershipRange(L,&Istart,&Iend);
115:   for (i=Istart;i<Iend;i++) {
116:     if (i>0) { MatSetValue(L,i,i-1,c,INSERT_VALUES); }
117:     if (i<n-1) { MatSetValue(L,i,i+1,c,INSERT_VALUES); }
118:     MatSetValue(L,i,i,-2.0*c,INSERT_VALUES);
119:   }
120:   MatAssemblyBegin(L,MAT_FINAL_ASSEMBLY);
121:   MatAssemblyEnd(L,MAT_FINAL_ASSEMBLY);
122:   MatCreateVecs(L,NULL,&u);
123:   VecDuplicate(u,&yex);
124:   VecDuplicate(u,&w);
125:   VecDuplicate(u,&z);

127:   /*
128:      Compute various vectors:
129:      - the exact solution yex = x*(1-x)*exp(tend)
130:      - the initial condition u = abs(x-0.5)-0.5
131:   */
132:   for (i=Istart;i<Iend;i++) {
133:     x = (i+1)*h;
134:     value = x*(1.0-x)*PetscExpReal(tend);
135:     VecSetValue(yex,i,value,INSERT_VALUES);
136:     value = PetscAbsReal(x-0.5)-0.5;
137:     VecSetValue(u,i,value,INSERT_VALUES);
138:   }
139:   VecAssemblyBegin(yex);
140:   VecAssemblyBegin(u);
141:   VecAssemblyEnd(yex);
142:   VecAssemblyEnd(u);
143:   VecViewFromOptions(yex,NULL,"-exact_sol");
144:   VecViewFromOptions(u,NULL,"-initial_cond");

146:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
147:               Create two MFN solvers, for exp() and phi_1()
148:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
149:   MFNCreate(PETSC_COMM_WORLD,&mfnexp);
150:   MFNSetOperator(mfnexp,L);
151:   MFNGetFN(mfnexp,&fexp);
152:   FNSetType(fexp,FNEXP);
153:   FNSetScale(fexp,deltat,1.0);
154:   MFNSetErrorIfNotConverged(mfnexp,PETSC_TRUE);
155:   MFNSetFromOptions(mfnexp);

157:   MFNCreate(PETSC_COMM_WORLD,&mfnphi);
158:   MFNSetOperator(mfnphi,L);
159:   MFNGetFN(mfnphi,&fphi);
160:   if (combine) {
161:     BuildFNPhi(fphi);
162:   } else {
163:     FNSetType(fphi,FNPHI);
164:     FNPhiSetIndex(fphi,1);
165:   }
166:   FNSetScale(fphi,deltat,1.0);
167:   MFNSetErrorIfNotConverged(mfnphi,PETSC_TRUE);
168:   MFNSetFromOptions(mfnphi);

170:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
171:              Solve the problem with the Norsett-Euler scheme
172:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
173:   steps = PetscRoundReal(tend/deltat);
174:   t = 0.0;
175:   for (k=0;k<steps;k++) {

177:     /* evaluate nonlinear part */
178:     VecGetArrayRead(u,&uarray);
179:     VecGetArray(w,&warray);
180:     for (i=Istart;i<Iend;i++) {
181:       x = (i+1)*h;
182:       uval = uarray[i-Istart];
183:       value = x*(1.0-x)*PetscExpReal(t);
184:       value = value + 2.0*PetscExpReal(t) - 1.0/(1.0+value*value);
185:       value = value + 1.0/(1.0+uval*uval);
186:       warray[i-Istart] = deltat*value;
187:     }
188:     VecRestoreArrayRead(u,&uarray);
189:     VecRestoreArray(w,&warray);
190:     MFNSolve(mfnphi,w,z);

192:     /* evaluate linear part */
193:     MFNSolve(mfnexp,u,u);
194:     VecAXPY(u,1.0,z);
195:     t = t + deltat;

197:   }
198:   VecViewFromOptions(u,NULL,"-computed_sol");

200:   /*
201:      Compare with exact solution and show error norm
202:   */
203:   VecCopy(u,z);
204:   VecAXPY(z,-1.0,yex);
205:   VecNorm(z,NORM_2,&nrmd);
206:   VecNorm(u,NORM_2,&nrmu);
207:   PetscPrintf(PETSC_COMM_WORLD," The relative error at t=%g is %.4f\n\n",(double)t,(double)(nrmd/nrmu));

209:   /*
210:      Free work space
211:   */
212:   MFNDestroy(&mfnexp);
213:   MFNDestroy(&mfnphi);
214:   MatDestroy(&L);
215:   VecDestroy(&u);
216:   VecDestroy(&yex);
217:   VecDestroy(&w);
218:   VecDestroy(&z);
219:   SlepcFinalize();
220:   return ierr;
221: }

223: /*TEST

225:    build:
226:       requires: c99

228:    test:
229:       suffix: 1
230:       args: -n 127 -tend 0.125 -mfn_tol 1e-3 -deltat 0.025
231:       timeoutfactor: 2

233:    test:
234:       suffix: 2
235:       args: -n 127 -tend 0.125 -mfn_tol 1e-3 -deltat 0.025 -combine
236:       filter: sed -e "s/ (combine)//"
237:       requires: !single
238:       output_file: output/ex39_1.out
239:       timeoutfactor: 2

241: TEST*/