Actual source code: ex2.c
1: /*
2: Formatted test for TS routines.
4: Solves U_t=F(t,u)
5: Where:
7: [2*u1+u2
8: F(t,u)= [u1+2*u2+u3
9: [ u2+2*u3
10: We can compare the solutions from euler, beuler and SUNDIALS to
11: see what is the difference.
13: */
15: static char help[] = "Solves a linear ODE. \n\n";
17: #include <petscts.h>
18: #include <petscpc.h>
20: extern PetscErrorCode RHSFunction(TS,PetscReal,Vec,Vec,void*);
21: extern PetscErrorCode RHSJacobian(TS,PetscReal,Vec,Mat,Mat,void*);
22: extern PetscErrorCode Monitor(TS,PetscInt,PetscReal,Vec,void*);
23: extern PetscErrorCode Initial(Vec,void*);
24: extern PetscErrorCode MyMatMult(Mat,Vec,Vec);
26: extern PetscReal solx(PetscReal);
27: extern PetscReal soly(PetscReal);
28: extern PetscReal solz(PetscReal);
30: int main(int argc,char **argv)
31: {
32: PetscInt time_steps = 100,steps;
33: Vec global;
34: PetscReal dt,ftime;
35: TS ts;
36: Mat A = 0,S;
38: PetscInitialize(&argc,&argv,(char*)0,help);
39: PetscOptionsGetInt(NULL,NULL,"-time",&time_steps,NULL);
41: /* set initial conditions */
42: VecCreate(PETSC_COMM_WORLD,&global);
43: VecSetSizes(global,PETSC_DECIDE,3);
44: VecSetFromOptions(global);
45: Initial(global,NULL);
47: /* make timestep context */
48: TSCreate(PETSC_COMM_WORLD,&ts);
49: TSSetProblemType(ts,TS_NONLINEAR);
50: TSMonitorSet(ts,Monitor,NULL,NULL);
51: dt = 0.001;
53: /*
54: The user provides the RHS and Jacobian
55: */
56: TSSetRHSFunction(ts,NULL,RHSFunction,NULL);
57: MatCreate(PETSC_COMM_WORLD,&A);
58: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,3,3);
59: MatSetFromOptions(A);
60: MatSetUp(A);
61: RHSJacobian(ts,0.0,global,A,A,NULL);
62: TSSetRHSJacobian(ts,A,A,RHSJacobian,NULL);
64: MatCreateShell(PETSC_COMM_WORLD,3,3,3,3,NULL,&S);
65: MatShellSetOperation(S,MATOP_MULT,(void (*)(void))MyMatMult);
66: TSSetRHSJacobian(ts,S,A,RHSJacobian,NULL);
68: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);
69: TSSetFromOptions(ts);
71: TSSetTimeStep(ts,dt);
72: TSSetMaxSteps(ts,time_steps);
73: TSSetMaxTime(ts,1);
74: TSSetSolution(ts,global);
76: TSSolve(ts,global);
77: TSGetSolveTime(ts,&ftime);
78: TSGetStepNumber(ts,&steps);
80: /* free the memories */
82: TSDestroy(&ts);
83: VecDestroy(&global);
84: MatDestroy(&A);
85: MatDestroy(&S);
87: PetscFinalize();
88: return 0;
89: }
91: PetscErrorCode MyMatMult(Mat S,Vec x,Vec y)
92: {
93: const PetscScalar *inptr;
94: PetscScalar *outptr;
96: VecGetArrayRead(x,&inptr);
97: VecGetArrayWrite(y,&outptr);
99: outptr[0] = 2.0*inptr[0]+inptr[1];
100: outptr[1] = inptr[0]+2.0*inptr[1]+inptr[2];
101: outptr[2] = inptr[1]+2.0*inptr[2];
103: VecRestoreArrayRead(x,&inptr);
104: VecRestoreArrayWrite(y,&outptr);
105: return 0;
106: }
108: /* this test problem has initial values (1,1,1). */
109: PetscErrorCode Initial(Vec global,void *ctx)
110: {
111: PetscScalar *localptr;
112: PetscInt i,mybase,myend,locsize;
114: /* determine starting point of each processor */
115: VecGetOwnershipRange(global,&mybase,&myend);
116: VecGetLocalSize(global,&locsize);
118: /* Initialize the array */
119: VecGetArrayWrite(global,&localptr);
120: for (i=0; i<locsize; i++) localptr[i] = 1.0;
122: if (mybase == 0) localptr[0]=1.0;
124: VecRestoreArrayWrite(global,&localptr);
125: return 0;
126: }
128: PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal time,Vec global,void *ctx)
129: {
130: VecScatter scatter;
131: IS from,to;
132: PetscInt i,n,*idx;
133: Vec tmp_vec;
134: PetscErrorCode ierr;
135: const PetscScalar *tmp;
137: /* Get the size of the vector */
138: VecGetSize(global,&n);
140: /* Set the index sets */
141: PetscMalloc1(n,&idx);
142: for (i=0; i<n; i++) idx[i]=i;
144: /* Create local sequential vectors */
145: VecCreateSeq(PETSC_COMM_SELF,n,&tmp_vec);
147: /* Create scatter context */
148: ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&from);
149: ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&to);
150: VecScatterCreate(global,from,tmp_vec,to,&scatter);
151: VecScatterBegin(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD);
152: VecScatterEnd(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD);
154: VecGetArrayRead(tmp_vec,&tmp);
155: PetscPrintf(PETSC_COMM_WORLD,"At t =%14.6e u = %14.6e %14.6e %14.6e \n",
156: (double)time,(double)PetscRealPart(tmp[0]),(double)PetscRealPart(tmp[1]),(double)PetscRealPart(tmp[2]));
157: PetscPrintf(PETSC_COMM_WORLD,"At t =%14.6e errors = %14.6e %14.6e %14.6e \n",
158: (double)time,(double)PetscRealPart(tmp[0]-solx(time)),(double)PetscRealPart(tmp[1]-soly(time)),(double)PetscRealPart(tmp[2]-solz(time)));
159: VecRestoreArrayRead(tmp_vec,&tmp);
160: VecScatterDestroy(&scatter);
161: ISDestroy(&from);
162: ISDestroy(&to);
163: PetscFree(idx);
164: VecDestroy(&tmp_vec);
165: return 0;
166: }
168: PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec globalin,Vec globalout,void *ctx)
169: {
170: PetscScalar *outptr;
171: const PetscScalar *inptr;
172: PetscInt i,n,*idx;
173: IS from,to;
174: VecScatter scatter;
175: Vec tmp_in,tmp_out;
177: /* Get the length of parallel vector */
178: VecGetSize(globalin,&n);
180: /* Set the index sets */
181: PetscMalloc1(n,&idx);
182: for (i=0; i<n; i++) idx[i]=i;
184: /* Create local sequential vectors */
185: VecCreateSeq(PETSC_COMM_SELF,n,&tmp_in);
186: VecDuplicate(tmp_in,&tmp_out);
188: /* Create scatter context */
189: ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&from);
190: ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&to);
191: VecScatterCreate(globalin,from,tmp_in,to,&scatter);
192: VecScatterBegin(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD);
193: VecScatterEnd(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD);
194: VecScatterDestroy(&scatter);
196: /*Extract income array */
197: VecGetArrayRead(tmp_in,&inptr);
199: /* Extract outcome array*/
200: VecGetArrayWrite(tmp_out,&outptr);
202: outptr[0] = 2.0*inptr[0]+inptr[1];
203: outptr[1] = inptr[0]+2.0*inptr[1]+inptr[2];
204: outptr[2] = inptr[1]+2.0*inptr[2];
206: VecRestoreArrayRead(tmp_in,&inptr);
207: VecRestoreArrayWrite(tmp_out,&outptr);
209: VecScatterCreate(tmp_out,from,globalout,to,&scatter);
210: VecScatterBegin(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD);
211: VecScatterEnd(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD);
213: /* Destroy idx aand scatter */
214: ISDestroy(&from);
215: ISDestroy(&to);
216: VecScatterDestroy(&scatter);
217: VecDestroy(&tmp_in);
218: VecDestroy(&tmp_out);
219: PetscFree(idx);
220: return 0;
221: }
223: PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec x,Mat A,Mat BB,void *ctx)
224: {
225: PetscScalar v[3];
226: const PetscScalar *tmp;
227: PetscInt idx[3],i;
229: idx[0]=0; idx[1]=1; idx[2]=2;
230: VecGetArrayRead(x,&tmp);
232: i = 0;
233: v[0] = 2.0; v[1] = 1.0; v[2] = 0.0;
234: MatSetValues(BB,1,&i,3,idx,v,INSERT_VALUES);
236: i = 1;
237: v[0] = 1.0; v[1] = 2.0; v[2] = 1.0;
238: MatSetValues(BB,1,&i,3,idx,v,INSERT_VALUES);
240: i = 2;
241: v[0] = 0.0; v[1] = 1.0; v[2] = 2.0;
242: MatSetValues(BB,1,&i,3,idx,v,INSERT_VALUES);
244: MatAssemblyBegin(BB,MAT_FINAL_ASSEMBLY);
245: MatAssemblyEnd(BB,MAT_FINAL_ASSEMBLY);
247: if (A != BB) {
248: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
249: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
250: }
251: VecRestoreArrayRead(x,&tmp);
253: return 0;
254: }
256: /*
257: The exact solutions
258: */
259: PetscReal solx(PetscReal t)
260: {
261: return PetscExpReal((2.0 - PetscSqrtReal(2.0))*t)/2.0 - PetscExpReal((2.0 - PetscSqrtReal(2.0))*t)/(2.0*PetscSqrtReal(2.0)) +
262: PetscExpReal((2.0 + PetscSqrtReal(2.0))*t)/2.0 + PetscExpReal((2.0 + PetscSqrtReal(2.0))*t)/(2.0*PetscSqrtReal(2.0));
263: }
265: PetscReal soly(PetscReal t)
266: {
267: return PetscExpReal((2.0 - PetscSqrtReal(2.0))*t)/2.0 - PetscExpReal((2.0 - PetscSqrtReal(2.0))*t)/PetscSqrtReal(2.0) +
268: PetscExpReal((2.0 + PetscSqrtReal(2.0))*t)/2.0 + PetscExpReal((2.0 + PetscSqrtReal(2.0))*t)/PetscSqrtReal(2.0);
269: }
271: PetscReal solz(PetscReal t)
272: {
273: return PetscExpReal((2.0 - PetscSqrtReal(2.0))*t)/2.0 - PetscExpReal((2.0 - PetscSqrtReal(2.0))*t)/(2.0*PetscSqrtReal(2.0)) +
274: PetscExpReal((2.0 + PetscSqrtReal(2.0))*t)/2.0 + PetscExpReal((2.0 + PetscSqrtReal(2.0))*t)/(2.0*PetscSqrtReal(2.0));
275: }
277: /*TEST
279: test:
280: suffix: euler
281: args: -ts_type euler
282: requires: !single
284: test:
285: suffix: beuler
286: args: -ts_type beuler
287: requires: !single
289: TEST*/