xref: /petsc/src/ts/tutorials/ex8.c (revision 030f984af8d8bb4c203755d35bded3c05b3d83ce)
1 
2 static char help[] = "Nonlinear DAE benchmark problems.\n";
3 
4 /*
5    Include "petscts.h" so that we can use TS solvers.  Note that this
6    file automatically includes:
7      petscsys.h       - base PETSc routines   petscvec.h - vectors
8      petscmat.h - matrices
9      petscis.h     - index sets            petscksp.h - Krylov subspace methods
10      petscviewer.h - viewers               petscpc.h  - preconditioners
11      petscksp.h   - linear solvers
12 */
13 #include <petscts.h>
14 
15 typedef struct _Problem* Problem;
16 struct _Problem {
17   PetscErrorCode (*destroy)(Problem);
18   TSIFunction    function;
19   TSIJacobian    jacobian;
20   PetscErrorCode (*solution)(PetscReal,Vec,void*);
21   MPI_Comm       comm;
22   PetscReal      final_time;
23   PetscInt       n;
24   PetscBool      hasexact;
25   void           *data;
26 };
27 
28 /*
29       Stiff 3-variable system from chemical reactions, due to Robertson (1966), problem ROBER in Hairer&Wanner, ODE 2, 1996
30 */
31 static PetscErrorCode RoberFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx)
32 {
33   PetscErrorCode    ierr;
34   PetscScalar       *f;
35   const PetscScalar *x,*xdot;
36 
37   PetscFunctionBeginUser;
38   ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr);
39   ierr = VecGetArrayRead(Xdot,&xdot);CHKERRQ(ierr);
40   ierr = VecGetArray(F,&f);CHKERRQ(ierr);
41   f[0] = xdot[0] + 0.04*x[0] - 1e4*x[1]*x[2];
42   f[1] = xdot[1] - 0.04*x[0] + 1e4*x[1]*x[2] + 3e7*PetscSqr(x[1]);
43   f[2] = xdot[2] - 3e7*PetscSqr(x[1]);
44   ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr);
45   ierr = VecRestoreArrayRead(Xdot,&xdot);CHKERRQ(ierr);
46   ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
47   PetscFunctionReturn(0);
48 }
49 
50 static PetscErrorCode RoberJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ctx)
51 {
52   PetscErrorCode    ierr;
53   PetscInt          rowcol[] = {0,1,2};
54   PetscScalar       J[3][3];
55   const PetscScalar *x,*xdot;
56 
57   PetscFunctionBeginUser;
58   ierr    = VecGetArrayRead(X,&x);CHKERRQ(ierr);
59   ierr    = VecGetArrayRead(Xdot,&xdot);CHKERRQ(ierr);
60   J[0][0] = a + 0.04;     J[0][1] = -1e4*x[2];                   J[0][2] = -1e4*x[1];
61   J[1][0] = -0.04;        J[1][1] = a + 1e4*x[2] + 3e7*2*x[1];   J[1][2] = 1e4*x[1];
62   J[2][0] = 0;            J[2][1] = -3e7*2*x[1];                 J[2][2] = a;
63   ierr    = MatSetValues(B,3,rowcol,3,rowcol,&J[0][0],INSERT_VALUES);CHKERRQ(ierr);
64   ierr    = VecRestoreArrayRead(X,&x);CHKERRQ(ierr);
65   ierr    = VecRestoreArrayRead(Xdot,&xdot);CHKERRQ(ierr);
66 
67   ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
68   ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
69   if (A != B) {
70     ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
71     ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
72   }
73   PetscFunctionReturn(0);
74 }
75 
76 static PetscErrorCode RoberSolution(PetscReal t,Vec X,void *ctx)
77 {
78   PetscErrorCode ierr;
79   PetscScalar    *x;
80 
81   PetscFunctionBeginUser;
82   if (t != 0) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"not implemented");
83   ierr = VecGetArray(X,&x);CHKERRQ(ierr);
84   x[0] = 1;
85   x[1] = 0;
86   x[2] = 0;
87   ierr = VecRestoreArray(X,&x);CHKERRQ(ierr);
88   PetscFunctionReturn(0);
89 }
90 
91 static PetscErrorCode RoberCreate(Problem p)
92 {
93 
94   PetscFunctionBeginUser;
95   p->destroy    = 0;
96   p->function   = &RoberFunction;
97   p->jacobian   = &RoberJacobian;
98   p->solution   = &RoberSolution;
99   p->final_time = 1e11;
100   p->n          = 3;
101   PetscFunctionReturn(0);
102 }
103 
104 /*
105      Stiff scalar valued problem
106 */
107 
108 typedef struct {
109   PetscReal lambda;
110 } CECtx;
111 
112 static PetscErrorCode CEDestroy(Problem p)
113 {
114   PetscErrorCode ierr;
115 
116   PetscFunctionBeginUser;
117   ierr = PetscFree(p->data);CHKERRQ(ierr);
118   PetscFunctionReturn(0);
119 }
120 
121 static PetscErrorCode CEFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx)
122 {
123   PetscErrorCode    ierr;
124   PetscReal         l = ((CECtx*)ctx)->lambda;
125   PetscScalar       *f;
126   const PetscScalar *x,*xdot;
127 
128   PetscFunctionBeginUser;
129   ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr);
130   ierr = VecGetArrayRead(Xdot,&xdot);CHKERRQ(ierr);
131   ierr = VecGetArray(F,&f);CHKERRQ(ierr);
132   f[0] = xdot[0] + l*(x[0] - PetscCosReal(t));
133 #if 0
134   ierr = PetscPrintf(PETSC_COMM_WORLD," f(t=%g,x=%g,xdot=%g) = %g\n",(double)t,(double)x[0],(double)xdot[0],(double)f[0]);CHKERRQ(ierr);
135 #endif
136   ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr);
137   ierr = VecRestoreArrayRead(Xdot,&xdot);CHKERRQ(ierr);
138   ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
139   PetscFunctionReturn(0);
140 }
141 
142 static PetscErrorCode CEJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ctx)
143 {
144   PetscReal         l = ((CECtx*)ctx)->lambda;
145   PetscErrorCode    ierr;
146   PetscInt          rowcol[] = {0};
147   PetscScalar       J[1][1];
148   const PetscScalar *x,*xdot;
149 
150   PetscFunctionBeginUser;
151   ierr    = VecGetArrayRead(X,&x);CHKERRQ(ierr);
152   ierr    = VecGetArrayRead(Xdot,&xdot);CHKERRQ(ierr);
153   J[0][0] = a + l;
154   ierr    = MatSetValues(B,1,rowcol,1,rowcol,&J[0][0],INSERT_VALUES);CHKERRQ(ierr);
155   ierr    = VecRestoreArrayRead(X,&x);CHKERRQ(ierr);
156   ierr    = VecRestoreArrayRead(Xdot,&xdot);CHKERRQ(ierr);
157 
158   ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
159   ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
160   if (A != B) {
161     ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
162     ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
163   }
164   PetscFunctionReturn(0);
165 }
166 
167 static PetscErrorCode CESolution(PetscReal t,Vec X,void *ctx)
168 {
169   PetscReal      l = ((CECtx*)ctx)->lambda;
170   PetscErrorCode ierr;
171   PetscScalar    *x;
172 
173   PetscFunctionBeginUser;
174   ierr = VecGetArray(X,&x);CHKERRQ(ierr);
175   x[0] = l/(l*l+1)*(l*PetscCosReal(t)+PetscSinReal(t)) - l*l/(l*l+1)*PetscExpReal(-l*t);
176   ierr = VecRestoreArray(X,&x);CHKERRQ(ierr);
177   PetscFunctionReturn(0);
178 }
179 
180 static PetscErrorCode CECreate(Problem p)
181 {
182   PetscErrorCode ierr;
183   CECtx          *ce;
184 
185   PetscFunctionBeginUser;
186   ierr    = PetscMalloc(sizeof(CECtx),&ce);CHKERRQ(ierr);
187   p->data = (void*)ce;
188 
189   p->destroy    = &CEDestroy;
190   p->function   = &CEFunction;
191   p->jacobian   = &CEJacobian;
192   p->solution   = &CESolution;
193   p->final_time = 10;
194   p->n          = 1;
195   p->hasexact   = PETSC_TRUE;
196 
197   ce->lambda = 10;
198   ierr       = PetscOptionsBegin(p->comm,NULL,"CE options","");CHKERRQ(ierr);
199   {
200     ierr = PetscOptionsReal("-problem_ce_lambda","Parameter controlling stiffness: xdot + lambda*(x - cos(t))","",ce->lambda,&ce->lambda,NULL);CHKERRQ(ierr);
201   }
202   ierr = PetscOptionsEnd();CHKERRQ(ierr);
203   PetscFunctionReturn(0);
204 }
205 
206 /*
207    Stiff 3-variable oscillatory system from chemical reactions. problem OREGO in Hairer&Wanner
208 */
209 static PetscErrorCode OregoFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx)
210 {
211   PetscErrorCode    ierr;
212   PetscScalar       *f;
213   const PetscScalar *x,*xdot;
214 
215   PetscFunctionBeginUser;
216   ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr);
217   ierr = VecGetArrayRead(Xdot,&xdot);CHKERRQ(ierr);
218   ierr = VecGetArray(F,&f);CHKERRQ(ierr);
219   f[0] = xdot[0] - 77.27*(x[1] + x[0]*(1. - 8.375e-6*x[0] - x[1]));
220   f[1] = xdot[1] - 1/77.27*(x[2] - (1. + x[0])*x[1]);
221   f[2] = xdot[2] - 0.161*(x[0] - x[2]);
222   ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr);
223   ierr = VecRestoreArrayRead(Xdot,&xdot);CHKERRQ(ierr);
224   ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
225   PetscFunctionReturn(0);
226 }
227 
228 static PetscErrorCode OregoJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ctx)
229 {
230   PetscErrorCode    ierr;
231   PetscInt          rowcol[] = {0,1,2};
232   PetscScalar       J[3][3];
233   const PetscScalar *x,*xdot;
234 
235   PetscFunctionBeginUser;
236   ierr    = VecGetArrayRead(X,&x);CHKERRQ(ierr);
237   ierr    = VecGetArrayRead(Xdot,&xdot);CHKERRQ(ierr);
238   J[0][0] = a - 77.27*((1. - 8.375e-6*x[0] - x[1]) - 8.375e-6*x[0]);
239   J[0][1] = -77.27*(1. - x[0]);
240   J[0][2] = 0;
241   J[1][0] = 1./77.27*x[1];
242   J[1][1] = a + 1./77.27*(1. + x[0]);
243   J[1][2] = -1./77.27;
244   J[2][0] = -0.161;
245   J[2][1] = 0;
246   J[2][2] = a + 0.161;
247   ierr    = MatSetValues(B,3,rowcol,3,rowcol,&J[0][0],INSERT_VALUES);CHKERRQ(ierr);
248   ierr    = VecRestoreArrayRead(X,&x);CHKERRQ(ierr);
249   ierr    = VecRestoreArrayRead(Xdot,&xdot);CHKERRQ(ierr);
250 
251   ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
252   ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
253   if (A != B) {
254     ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
255     ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
256   }
257   PetscFunctionReturn(0);
258 }
259 
260 static PetscErrorCode OregoSolution(PetscReal t,Vec X,void *ctx)
261 {
262   PetscErrorCode ierr;
263   PetscScalar    *x;
264 
265   PetscFunctionBeginUser;
266   if (t != 0) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"not implemented");
267   ierr = VecGetArray(X,&x);CHKERRQ(ierr);
268   x[0] = 1;
269   x[1] = 2;
270   x[2] = 3;
271   ierr = VecRestoreArray(X,&x);CHKERRQ(ierr);
272   PetscFunctionReturn(0);
273 }
274 
275 static PetscErrorCode OregoCreate(Problem p)
276 {
277 
278   PetscFunctionBeginUser;
279   p->destroy    = 0;
280   p->function   = &OregoFunction;
281   p->jacobian   = &OregoJacobian;
282   p->solution   = &OregoSolution;
283   p->final_time = 360;
284   p->n          = 3;
285   PetscFunctionReturn(0);
286 }
287 
288 /*
289    User-defined monitor for comparing to exact solutions when possible
290 */
291 typedef struct {
292   MPI_Comm comm;
293   Problem  problem;
294   Vec      x;
295 } MonitorCtx;
296 
297 static PetscErrorCode MonitorError(TS ts,PetscInt step,PetscReal t,Vec x,void *ctx)
298 {
299   PetscErrorCode ierr;
300   MonitorCtx     *mon = (MonitorCtx*)ctx;
301   PetscReal      h,nrm_x,nrm_exact,nrm_diff;
302 
303   PetscFunctionBeginUser;
304   if (!mon->problem->solution) PetscFunctionReturn(0);
305   ierr = (*mon->problem->solution)(t,mon->x,mon->problem->data);CHKERRQ(ierr);
306   ierr = VecNorm(x,NORM_2,&nrm_x);CHKERRQ(ierr);
307   ierr = VecNorm(mon->x,NORM_2,&nrm_exact);CHKERRQ(ierr);
308   ierr = VecAYPX(mon->x,-1,x);CHKERRQ(ierr);
309   ierr = VecNorm(mon->x,NORM_2,&nrm_diff);CHKERRQ(ierr);
310   ierr = TSGetTimeStep(ts,&h);CHKERRQ(ierr);
311   if (step < 0) {
312     ierr = PetscPrintf(mon->comm,"Interpolated final solution ");CHKERRQ(ierr);
313   }
314   ierr = PetscPrintf(mon->comm,"step %4D t=%12.8e h=% 8.2e  |x|=%9.2e  |x_e|=%9.2e  |x-x_e|=%9.2e\n",step,(double)t,(double)h,(double)nrm_x,(double)nrm_exact,(double)nrm_diff);CHKERRQ(ierr);
315   PetscFunctionReturn(0);
316 }
317 
318 int main(int argc,char **argv)
319 {
320   PetscFunctionList plist = NULL;
321   char              pname[256];
322   TS                ts;            /* nonlinear solver */
323   Vec               x,r;           /* solution, residual vectors */
324   Mat               A;             /* Jacobian matrix */
325   Problem           problem;
326   PetscBool         use_monitor = PETSC_FALSE;
327   PetscBool         use_result = PETSC_FALSE;
328   PetscInt          steps,nonlinits,linits,snesfails,rejects;
329   PetscReal         ftime;
330   MonitorCtx        mon;
331   PetscErrorCode    ierr;
332   PetscMPIInt       size;
333 
334   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
335      Initialize program
336      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
337   ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
338   ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRMPI(ierr);
339   if (size > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"Only for sequential runs");
340 
341   /* Register the available problems */
342   ierr = PetscFunctionListAdd(&plist,"rober",&RoberCreate);CHKERRQ(ierr);
343   ierr = PetscFunctionListAdd(&plist,"ce",&CECreate);CHKERRQ(ierr);
344   ierr = PetscFunctionListAdd(&plist,"orego",&OregoCreate);CHKERRQ(ierr);
345   ierr = PetscStrcpy(pname,"ce");CHKERRQ(ierr);
346 
347   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
348     Set runtime options
349     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
350   ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Timestepping benchmark options","");CHKERRQ(ierr);
351   {
352     ierr        = PetscOptionsFList("-problem_type","Name of problem to run","",plist,pname,pname,sizeof(pname),NULL);CHKERRQ(ierr);
353     use_monitor = PETSC_FALSE;
354     ierr        = PetscOptionsBool("-monitor_error","Display errors relative to exact solutions","",use_monitor,&use_monitor,NULL);CHKERRQ(ierr);
355     ierr        = PetscOptionsBool("-monitor_result","Display result","",use_result,&use_result,NULL);CHKERRQ(ierr);
356   }
357   ierr = PetscOptionsEnd();CHKERRQ(ierr);
358 
359   /* Create the new problem */
360   ierr          = PetscNew(&problem);CHKERRQ(ierr);
361   problem->comm = MPI_COMM_WORLD;
362   {
363     PetscErrorCode (*pcreate)(Problem);
364 
365     ierr = PetscFunctionListFind(plist,pname,&pcreate);CHKERRQ(ierr);
366     if (!pcreate) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"No problem '%s'",pname);
367     ierr = (*pcreate)(problem);CHKERRQ(ierr);
368   }
369 
370   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
371     Create necessary matrix and vectors
372     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
373   ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
374   ierr = MatSetSizes(A,problem->n,problem->n,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr);
375   ierr = MatSetFromOptions(A);CHKERRQ(ierr);
376   ierr = MatSetUp(A);CHKERRQ(ierr);
377 
378   ierr = MatCreateVecs(A,&x,NULL);CHKERRQ(ierr);
379   ierr = VecDuplicate(x,&r);CHKERRQ(ierr);
380 
381   mon.comm    = PETSC_COMM_WORLD;
382   mon.problem = problem;
383   ierr        = VecDuplicate(x,&mon.x);CHKERRQ(ierr);
384 
385   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
386      Create timestepping solver context
387      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
388   ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr);
389   ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr);
390   ierr = TSSetType(ts,TSROSW);CHKERRQ(ierr); /* Rosenbrock-W */
391   ierr = TSSetIFunction(ts,NULL,problem->function,problem->data);CHKERRQ(ierr);
392   ierr = TSSetIJacobian(ts,A,A,problem->jacobian,problem->data);CHKERRQ(ierr);
393   ierr = TSSetMaxTime(ts,problem->final_time);CHKERRQ(ierr);
394   ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr);
395   ierr = TSSetMaxStepRejections(ts,10);CHKERRQ(ierr);
396   ierr = TSSetMaxSNESFailures(ts,-1);CHKERRQ(ierr); /* unlimited */
397   if (use_monitor) {
398     ierr = TSMonitorSet(ts,&MonitorError,&mon,NULL);CHKERRQ(ierr);
399   }
400 
401   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
402      Set initial conditions
403    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
404   ierr = (*problem->solution)(0,x,problem->data);CHKERRQ(ierr);
405   ierr = TSSetTimeStep(ts,.001);CHKERRQ(ierr);
406   ierr = TSSetSolution(ts,x);CHKERRQ(ierr);
407 
408   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
409      Set runtime options
410    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
411   ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
412 
413   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
414      Solve nonlinear system
415      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
416   ierr = TSSolve(ts,x);CHKERRQ(ierr);
417   ierr = TSGetSolveTime(ts,&ftime);CHKERRQ(ierr);
418   ierr = TSGetStepNumber(ts,&steps);CHKERRQ(ierr);
419   ierr = TSGetSNESFailures(ts,&snesfails);CHKERRQ(ierr);
420   ierr = TSGetStepRejections(ts,&rejects);CHKERRQ(ierr);
421   ierr = TSGetSNESIterations(ts,&nonlinits);CHKERRQ(ierr);
422   ierr = TSGetKSPIterations(ts,&linits);CHKERRQ(ierr);
423   if (use_result) {
424     ierr = PetscPrintf(PETSC_COMM_WORLD,"steps %D (%D rejected, %D SNES fails), ftime %g, nonlinits %D, linits %D\n",steps,rejects,snesfails,(double)ftime,nonlinits,linits);CHKERRQ(ierr);
425   }
426 
427   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
428      Free work space.  All PETSc objects should be destroyed when they
429      are no longer needed.
430    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
431   ierr = MatDestroy(&A);CHKERRQ(ierr);
432   ierr = VecDestroy(&x);CHKERRQ(ierr);
433   ierr = VecDestroy(&r);CHKERRQ(ierr);
434   ierr = VecDestroy(&mon.x);CHKERRQ(ierr);
435   ierr = TSDestroy(&ts);CHKERRQ(ierr);
436   if (problem->destroy) {
437     ierr = (*problem->destroy)(problem);CHKERRQ(ierr);
438   }
439   ierr = PetscFree(problem);CHKERRQ(ierr);
440   ierr = PetscFunctionListDestroy(&plist);CHKERRQ(ierr);
441 
442   ierr = PetscFinalize();
443   return ierr;
444 }
445 
446 /*TEST
447 
448     test:
449       requires: !complex
450       args:  -monitor_result -monitor_error -ts_atol 1e-2 -ts_rtol 1e-2 -ts_exact_final_time interpolate -ts_type arkimex
451 
452     test:
453       suffix: 2
454       requires: !single !complex
455       args: -monitor_result -ts_atol 1e-2 -ts_rtol 1e-2 -ts_max_time 15 -ts_type arkimex -ts_arkimex_type 2e -problem_type orego -ts_arkimex_initial_guess_extrapolate 0 -ts_adapt_time_step_increase_delay 4
456 
457     test:
458       suffix: 3
459       requires: !single !complex
460       args: -monitor_result -ts_atol 1e-2 -ts_rtol 1e-2 -ts_max_time 15 -ts_type arkimex -ts_arkimex_type 2e -problem_type orego -ts_arkimex_initial_guess_extrapolate 1
461 
462     test:
463       suffix: 4
464 
465     test:
466       suffix: 5
467       args: -snes_lag_jacobian 20 -snes_lag_jacobian_persists
468 
469 TEST*/
470