Actual source code: ex2.c
petsc-3.11.4 2019-09-28
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: {
33: PetscInt time_steps = 100,steps;
34: PetscMPIInt size;
35: Vec global;
36: PetscReal dt,ftime;
37: TS ts;
38: Mat A = 0,S;
40: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
41: MPI_Comm_size(PETSC_COMM_WORLD,&size);
43: PetscOptionsGetInt(NULL,NULL,"-time",&time_steps,NULL);
45: /* set initial conditions */
46: VecCreate(PETSC_COMM_WORLD,&global);
47: VecSetSizes(global,PETSC_DECIDE,3);
48: VecSetFromOptions(global);
49: Initial(global,NULL);
51: /* make timestep context */
52: TSCreate(PETSC_COMM_WORLD,&ts);
53: TSSetProblemType(ts,TS_NONLINEAR);
54: TSMonitorSet(ts,Monitor,NULL,NULL);
56: dt = 0.001;
58: /*
59: The user provides the RHS and Jacobian
60: */
61: TSSetRHSFunction(ts,NULL,RHSFunction,NULL);
62: MatCreate(PETSC_COMM_WORLD,&A);
63: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,3,3);
64: MatSetFromOptions(A);
65: MatSetUp(A);
66: RHSJacobian(ts,0.0,global,A,A,NULL);
67: TSSetRHSJacobian(ts,A,A,RHSJacobian,NULL);
69: MatCreateShell(PETSC_COMM_WORLD,3,3,3,3,NULL,&S);
70: MatShellSetOperation(S,MATOP_MULT,(void (*)(void))MyMatMult);
71: TSSetRHSJacobian(ts,S,A,RHSJacobian,NULL);
73: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);
74: TSSetFromOptions(ts);
76: TSSetTimeStep(ts,dt);
77: TSSetMaxSteps(ts,time_steps);
78: TSSetMaxTime(ts,1);
79: TSSetSolution(ts,global);
81: TSSolve(ts,global);
82: TSGetSolveTime(ts,&ftime);
83: TSGetStepNumber(ts,&steps);
86: /* free the memories */
88: TSDestroy(&ts);
89: VecDestroy(&global);
90: MatDestroy(&A);
91: MatDestroy(&S);
93: PetscFinalize();
94: return ierr;
95: }
97: PetscErrorCode MyMatMult(Mat S,Vec x,Vec y)
98: {
99: PetscErrorCode ierr;
100: const PetscScalar *inptr;
101: PetscScalar *outptr;
104: VecGetArrayRead(x,&inptr);
105: VecGetArray(y,&outptr);
107: outptr[0] = 2.0*inptr[0]+inptr[1];
108: outptr[1] = inptr[0]+2.0*inptr[1]+inptr[2];
109: outptr[2] = inptr[1]+2.0*inptr[2];
111: VecRestoreArrayRead(x,&inptr);
112: VecRestoreArray(y,&outptr);
113: return(0);
114: }
117: /* -------------------------------------------------------------------*/
118: /* this test problem has initial values (1,1,1). */
119: PetscErrorCode Initial(Vec global,void *ctx)
120: {
121: PetscScalar *localptr;
122: PetscInt i,mybase,myend,locsize;
125: /* determine starting point of each processor */
126: VecGetOwnershipRange(global,&mybase,&myend);
127: VecGetLocalSize(global,&locsize);
129: /* Initialize the array */
130: VecGetArray(global,&localptr);
131: for (i=0; i<locsize; i++) localptr[i] = 1.0;
133: if (mybase == 0) localptr[0]=1.0;
135: VecRestoreArray(global,&localptr);
136: return 0;
137: }
139: PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal time,Vec global,void *ctx)
140: {
141: VecScatter scatter;
142: IS from,to;
143: PetscInt i,n,*idx;
144: Vec tmp_vec;
145: PetscErrorCode ierr;
146: const PetscScalar *tmp;
148: /* Get the size of the vector */
149: VecGetSize(global,&n);
151: /* Set the index sets */
152: PetscMalloc1(n,&idx);
153: for (i=0; i<n; i++) idx[i]=i;
155: /* Create local sequential vectors */
156: VecCreateSeq(PETSC_COMM_SELF,n,&tmp_vec);
158: /* Create scatter context */
159: ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&from);
160: ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&to);
161: VecScatterCreate(global,from,tmp_vec,to,&scatter);
162: VecScatterBegin(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD);
163: VecScatterEnd(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD);
165: VecGetArrayRead(tmp_vec,&tmp);
166: PetscPrintf(PETSC_COMM_WORLD,"At t =%14.6e u = %14.6e %14.6e %14.6e \n",
167: (double)time,(double)PetscRealPart(tmp[0]),(double)PetscRealPart(tmp[1]),(double)PetscRealPart(tmp[2]));
168: PetscPrintf(PETSC_COMM_WORLD,"At t =%14.6e errors = %14.6e %14.6e %14.6e \n",
169: (double)time,(double)PetscRealPart(tmp[0]-solx(time)),(double)PetscRealPart(tmp[1]-soly(time)),(double)PetscRealPart(tmp[2]-solz(time)));
170: VecRestoreArrayRead(tmp_vec,&tmp);
171: VecScatterDestroy(&scatter);
172: ISDestroy(&from);
173: ISDestroy(&to);
174: PetscFree(idx);
175: VecDestroy(&tmp_vec);
176: return 0;
177: }
179: PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec globalin,Vec globalout,void *ctx)
180: {
181: PetscScalar *outptr;
182: const PetscScalar *inptr;
183: PetscInt i,n,*idx;
184: PetscErrorCode ierr;
185: IS from,to;
186: VecScatter scatter;
187: Vec tmp_in,tmp_out;
189: /* Get the length of parallel vector */
190: VecGetSize(globalin,&n);
192: /* Set the index sets */
193: PetscMalloc1(n,&idx);
194: for (i=0; i<n; i++) idx[i]=i;
196: /* Create local sequential vectors */
197: VecCreateSeq(PETSC_COMM_SELF,n,&tmp_in);
198: VecDuplicate(tmp_in,&tmp_out);
200: /* Create scatter context */
201: ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&from);
202: ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&to);
203: VecScatterCreate(globalin,from,tmp_in,to,&scatter);
204: VecScatterBegin(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD);
205: VecScatterEnd(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD);
206: VecScatterDestroy(&scatter);
208: /*Extract income array */
209: VecGetArrayRead(tmp_in,&inptr);
211: /* Extract outcome array*/
212: VecGetArray(tmp_out,&outptr);
214: outptr[0] = 2.0*inptr[0]+inptr[1];
215: outptr[1] = inptr[0]+2.0*inptr[1]+inptr[2];
216: outptr[2] = inptr[1]+2.0*inptr[2];
218: VecRestoreArrayRead(tmp_in,&inptr);
219: VecRestoreArray(tmp_out,&outptr);
221: VecScatterCreate(tmp_out,from,globalout,to,&scatter);
222: VecScatterBegin(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD);
223: VecScatterEnd(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD);
225: /* Destroy idx aand scatter */
226: ISDestroy(&from);
227: ISDestroy(&to);
228: VecScatterDestroy(&scatter);
229: VecDestroy(&tmp_in);
230: VecDestroy(&tmp_out);
231: PetscFree(idx);
232: return 0;
233: }
235: PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec x,Mat A,Mat BB,void *ctx)
236: {
237: PetscScalar v[3];
238: const PetscScalar *tmp;
239: PetscInt idx[3],i;
240: PetscErrorCode ierr;
242: idx[0]=0; idx[1]=1; idx[2]=2;
243: VecGetArrayRead(x,&tmp);
245: i = 0;
246: v[0] = 2.0; v[1] = 1.0; v[2] = 0.0;
247: MatSetValues(BB,1,&i,3,idx,v,INSERT_VALUES);
249: i = 1;
250: v[0] = 1.0; v[1] = 2.0; v[2] = 1.0;
251: MatSetValues(BB,1,&i,3,idx,v,INSERT_VALUES);
253: i = 2;
254: v[0] = 0.0; v[1] = 1.0; v[2] = 2.0;
255: MatSetValues(BB,1,&i,3,idx,v,INSERT_VALUES);
257: MatAssemblyBegin(BB,MAT_FINAL_ASSEMBLY);
258: MatAssemblyEnd(BB,MAT_FINAL_ASSEMBLY);
260: if (A != BB) {
261: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
262: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
263: }
264: VecRestoreArrayRead(x,&tmp);
266: return 0;
267: }
269: /*
270: The exact solutions
271: */
272: PetscReal solx(PetscReal t)
273: {
274: return PetscExpReal((2.0 - PetscSqrtReal(2.0))*t)/2.0 - PetscExpReal((2.0 - PetscSqrtReal(2.0))*t)/(2.0*PetscSqrtReal(2.0)) +
275: PetscExpReal((2.0 + PetscSqrtReal(2.0))*t)/2.0 + PetscExpReal((2.0 + PetscSqrtReal(2.0))*t)/(2.0*PetscSqrtReal(2.0));
276: }
278: PetscReal soly(PetscReal t)
279: {
280: return PetscExpReal((2.0 - PetscSqrtReal(2.0))*t)/2.0 - PetscExpReal((2.0 - PetscSqrtReal(2.0))*t)/PetscSqrtReal(2.0) +
281: PetscExpReal((2.0 + PetscSqrtReal(2.0))*t)/2.0 + PetscExpReal((2.0 + PetscSqrtReal(2.0))*t)/PetscSqrtReal(2.0);
282: }
284: PetscReal solz(PetscReal t)
285: {
286: return PetscExpReal((2.0 - PetscSqrtReal(2.0))*t)/2.0 - PetscExpReal((2.0 - PetscSqrtReal(2.0))*t)/(2.0*PetscSqrtReal(2.0)) +
287: PetscExpReal((2.0 + PetscSqrtReal(2.0))*t)/2.0 + PetscExpReal((2.0 + PetscSqrtReal(2.0))*t)/(2.0*PetscSqrtReal(2.0));
288: }
291: /*TEST
293: test:
294: suffix: euler
295: args: -ts_type euler
296: requires: !single
298: test:
299: suffix: beuler
300: args: -ts_type beuler
301: requires: !single
303: TEST*/