Actual source code: ex2.c
petsc-3.14.6 2021-03-30
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: Vec global;
35: PetscReal dt,ftime;
36: TS ts;
37: Mat A = 0,S;
39: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
40: PetscOptionsGetInt(NULL,NULL,"-time",&time_steps,NULL);
42: /* set initial conditions */
43: VecCreate(PETSC_COMM_WORLD,&global);
44: VecSetSizes(global,PETSC_DECIDE,3);
45: VecSetFromOptions(global);
46: Initial(global,NULL);
48: /* make timestep context */
49: TSCreate(PETSC_COMM_WORLD,&ts);
50: TSSetProblemType(ts,TS_NONLINEAR);
51: TSMonitorSet(ts,Monitor,NULL,NULL);
52: dt = 0.001;
54: /*
55: The user provides the RHS and Jacobian
56: */
57: TSSetRHSFunction(ts,NULL,RHSFunction,NULL);
58: MatCreate(PETSC_COMM_WORLD,&A);
59: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,3,3);
60: MatSetFromOptions(A);
61: MatSetUp(A);
62: RHSJacobian(ts,0.0,global,A,A,NULL);
63: TSSetRHSJacobian(ts,A,A,RHSJacobian,NULL);
65: MatCreateShell(PETSC_COMM_WORLD,3,3,3,3,NULL,&S);
66: MatShellSetOperation(S,MATOP_MULT,(void (*)(void))MyMatMult);
67: TSSetRHSJacobian(ts,S,A,RHSJacobian,NULL);
69: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);
70: TSSetFromOptions(ts);
72: TSSetTimeStep(ts,dt);
73: TSSetMaxSteps(ts,time_steps);
74: TSSetMaxTime(ts,1);
75: TSSetSolution(ts,global);
77: TSSolve(ts,global);
78: TSGetSolveTime(ts,&ftime);
79: TSGetStepNumber(ts,&steps);
81: /* free the memories */
83: TSDestroy(&ts);
84: VecDestroy(&global);
85: MatDestroy(&A);
86: MatDestroy(&S);
88: PetscFinalize();
89: return ierr;
90: }
92: PetscErrorCode MyMatMult(Mat S,Vec x,Vec y)
93: {
94: PetscErrorCode ierr;
95: const PetscScalar *inptr;
96: PetscScalar *outptr;
99: VecGetArrayRead(x,&inptr);
100: VecGetArrayWrite(y,&outptr);
102: outptr[0] = 2.0*inptr[0]+inptr[1];
103: outptr[1] = inptr[0]+2.0*inptr[1]+inptr[2];
104: outptr[2] = inptr[1]+2.0*inptr[2];
106: VecRestoreArrayRead(x,&inptr);
107: VecRestoreArrayWrite(y,&outptr);
108: return(0);
109: }
111: /* this test problem has initial values (1,1,1). */
112: PetscErrorCode Initial(Vec global,void *ctx)
113: {
114: PetscScalar *localptr;
115: PetscInt i,mybase,myend,locsize;
118: /* determine starting point of each processor */
119: VecGetOwnershipRange(global,&mybase,&myend);
120: VecGetLocalSize(global,&locsize);
122: /* Initialize the array */
123: VecGetArrayWrite(global,&localptr);
124: for (i=0; i<locsize; i++) localptr[i] = 1.0;
126: if (mybase == 0) localptr[0]=1.0;
128: VecRestoreArrayWrite(global,&localptr);
129: return 0;
130: }
132: PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal time,Vec global,void *ctx)
133: {
134: VecScatter scatter;
135: IS from,to;
136: PetscInt i,n,*idx;
137: Vec tmp_vec;
138: PetscErrorCode ierr;
139: const PetscScalar *tmp;
141: /* Get the size of the vector */
142: VecGetSize(global,&n);
144: /* Set the index sets */
145: PetscMalloc1(n,&idx);
146: for (i=0; i<n; i++) idx[i]=i;
148: /* Create local sequential vectors */
149: VecCreateSeq(PETSC_COMM_SELF,n,&tmp_vec);
151: /* Create scatter context */
152: ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&from);
153: ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&to);
154: VecScatterCreate(global,from,tmp_vec,to,&scatter);
155: VecScatterBegin(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD);
156: VecScatterEnd(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD);
158: VecGetArrayRead(tmp_vec,&tmp);
159: PetscPrintf(PETSC_COMM_WORLD,"At t =%14.6e u = %14.6e %14.6e %14.6e \n",
160: (double)time,(double)PetscRealPart(tmp[0]),(double)PetscRealPart(tmp[1]),(double)PetscRealPart(tmp[2]));
161: PetscPrintf(PETSC_COMM_WORLD,"At t =%14.6e errors = %14.6e %14.6e %14.6e \n",
162: (double)time,(double)PetscRealPart(tmp[0]-solx(time)),(double)PetscRealPart(tmp[1]-soly(time)),(double)PetscRealPart(tmp[2]-solz(time)));
163: VecRestoreArrayRead(tmp_vec,&tmp);
164: VecScatterDestroy(&scatter);
165: ISDestroy(&from);
166: ISDestroy(&to);
167: PetscFree(idx);
168: VecDestroy(&tmp_vec);
169: return 0;
170: }
172: PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec globalin,Vec globalout,void *ctx)
173: {
174: PetscScalar *outptr;
175: const PetscScalar *inptr;
176: PetscInt i,n,*idx;
177: PetscErrorCode ierr;
178: IS from,to;
179: VecScatter scatter;
180: Vec tmp_in,tmp_out;
182: /* Get the length of parallel vector */
183: VecGetSize(globalin,&n);
185: /* Set the index sets */
186: PetscMalloc1(n,&idx);
187: for (i=0; i<n; i++) idx[i]=i;
189: /* Create local sequential vectors */
190: VecCreateSeq(PETSC_COMM_SELF,n,&tmp_in);
191: VecDuplicate(tmp_in,&tmp_out);
193: /* Create scatter context */
194: ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&from);
195: ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&to);
196: VecScatterCreate(globalin,from,tmp_in,to,&scatter);
197: VecScatterBegin(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD);
198: VecScatterEnd(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD);
199: VecScatterDestroy(&scatter);
201: /*Extract income array */
202: VecGetArrayRead(tmp_in,&inptr);
204: /* Extract outcome array*/
205: VecGetArrayWrite(tmp_out,&outptr);
207: outptr[0] = 2.0*inptr[0]+inptr[1];
208: outptr[1] = inptr[0]+2.0*inptr[1]+inptr[2];
209: outptr[2] = inptr[1]+2.0*inptr[2];
211: VecRestoreArrayRead(tmp_in,&inptr);
212: VecRestoreArrayWrite(tmp_out,&outptr);
214: VecScatterCreate(tmp_out,from,globalout,to,&scatter);
215: VecScatterBegin(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD);
216: VecScatterEnd(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD);
218: /* Destroy idx aand scatter */
219: ISDestroy(&from);
220: ISDestroy(&to);
221: VecScatterDestroy(&scatter);
222: VecDestroy(&tmp_in);
223: VecDestroy(&tmp_out);
224: PetscFree(idx);
225: return 0;
226: }
228: PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec x,Mat A,Mat BB,void *ctx)
229: {
230: PetscScalar v[3];
231: const PetscScalar *tmp;
232: PetscInt idx[3],i;
233: PetscErrorCode ierr;
235: idx[0]=0; idx[1]=1; idx[2]=2;
236: VecGetArrayRead(x,&tmp);
238: i = 0;
239: v[0] = 2.0; v[1] = 1.0; v[2] = 0.0;
240: MatSetValues(BB,1,&i,3,idx,v,INSERT_VALUES);
242: i = 1;
243: v[0] = 1.0; v[1] = 2.0; v[2] = 1.0;
244: MatSetValues(BB,1,&i,3,idx,v,INSERT_VALUES);
246: i = 2;
247: v[0] = 0.0; v[1] = 1.0; v[2] = 2.0;
248: MatSetValues(BB,1,&i,3,idx,v,INSERT_VALUES);
250: MatAssemblyBegin(BB,MAT_FINAL_ASSEMBLY);
251: MatAssemblyEnd(BB,MAT_FINAL_ASSEMBLY);
253: if (A != BB) {
254: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
255: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
256: }
257: VecRestoreArrayRead(x,&tmp);
259: return 0;
260: }
262: /*
263: The exact solutions
264: */
265: PetscReal solx(PetscReal t)
266: {
267: return PetscExpReal((2.0 - PetscSqrtReal(2.0))*t)/2.0 - PetscExpReal((2.0 - PetscSqrtReal(2.0))*t)/(2.0*PetscSqrtReal(2.0)) +
268: PetscExpReal((2.0 + PetscSqrtReal(2.0))*t)/2.0 + PetscExpReal((2.0 + PetscSqrtReal(2.0))*t)/(2.0*PetscSqrtReal(2.0));
269: }
271: PetscReal soly(PetscReal t)
272: {
273: return PetscExpReal((2.0 - PetscSqrtReal(2.0))*t)/2.0 - PetscExpReal((2.0 - PetscSqrtReal(2.0))*t)/PetscSqrtReal(2.0) +
274: PetscExpReal((2.0 + PetscSqrtReal(2.0))*t)/2.0 + PetscExpReal((2.0 + PetscSqrtReal(2.0))*t)/PetscSqrtReal(2.0);
275: }
277: PetscReal solz(PetscReal t)
278: {
279: return PetscExpReal((2.0 - PetscSqrtReal(2.0))*t)/2.0 - PetscExpReal((2.0 - PetscSqrtReal(2.0))*t)/(2.0*PetscSqrtReal(2.0)) +
280: PetscExpReal((2.0 + PetscSqrtReal(2.0))*t)/2.0 + PetscExpReal((2.0 + PetscSqrtReal(2.0))*t)/(2.0*PetscSqrtReal(2.0));
281: }
284: /*TEST
286: test:
287: suffix: euler
288: args: -ts_type euler
289: requires: !single
291: test:
292: suffix: beuler
293: args: -ts_type beuler
294: requires: !single
296: TEST*/