Actual source code: ex2.c
1: /*
2: Formatted test for TS routines.
4: Solves U_t=F(t,u)
5: Where:
6:
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 nonlinear ODE. \n\n";
17: #include petscsys.h
18: #include petscts.h
19: #include petscpc.h
32: int main(int argc,char **argv)
33: {
35: PetscInt time_steps = 100,steps;
36: PetscMPIInt size;
37: Vec global;
38: PetscReal dt,ftime;
39: TS ts;
40: MatStructure A_structure;
41: Mat A = 0;
42:
43: PetscInitialize(&argc,&argv,(char*)0,help);
44: MPI_Comm_size(PETSC_COMM_WORLD,&size);
45:
46: PetscOptionsGetInt(PETSC_NULL,"-time",&time_steps,PETSC_NULL);
47:
48: /* set initial conditions */
49: VecCreate(PETSC_COMM_WORLD,&global);
50: VecSetSizes(global,PETSC_DECIDE,3);
51: VecSetFromOptions(global);
52: Initial(global,PETSC_NULL);
53:
54: /* make timestep context */
55: TSCreate(PETSC_COMM_WORLD,&ts);
56: TSSetProblemType(ts,TS_NONLINEAR);
57: TSMonitorSet(ts,Monitor,PETSC_NULL,PETSC_NULL);
59: dt = 0.1;
61: /*
62: The user provides the RHS and Jacobian
63: */
64: TSSetRHSFunction(ts,RHSFunction,NULL);
65: MatCreate(PETSC_COMM_WORLD,&A);
66: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,3,3);
67: MatSetFromOptions(A);
68: RHSJacobian(ts,0.0,global,&A,&A,&A_structure,NULL);
69: TSSetRHSJacobian(ts,A,A,RHSJacobian,NULL);
70:
71: TSSetFromOptions(ts);
73: TSSetInitialTimeStep(ts,0.0,dt);
74: TSSetDuration(ts,time_steps,1);
75: TSSetSolution(ts,global);
78: TSSetUp(ts);
79: TSStep(ts,&steps,&ftime);
82: /* free the memories */
83:
84: TSDestroy(ts);
85: VecDestroy(global);
86: ierr= MatDestroy(A);
88: PetscFinalize();
89: return 0;
90: }
92: /* -------------------------------------------------------------------*/
95: /* this test problem has initial values (1,1,1). */
96: PetscErrorCode Initial(Vec global,void *ctx)
97: {
98: PetscScalar *localptr;
99: PetscInt i,mybase,myend,locsize;
102: /* determine starting point of each processor */
103: VecGetOwnershipRange(global,&mybase,&myend);
104: VecGetLocalSize(global,&locsize);
106: /* Initialize the array */
107: VecGetArray(global,&localptr);
108: for (i=0; i<locsize; i++) {
109: localptr[i] = 1.0;
110: }
111:
112: if (mybase == 0) localptr[0]=1.0;
114: VecRestoreArray(global,&localptr);
115: return 0;
116: }
120: PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal time,Vec global,void *ctx)
121: {
122: VecScatter scatter;
123: IS from,to;
124: PetscInt i,n,*idx;
125: Vec tmp_vec;
127: PetscScalar *tmp;
129: /* Get the size of the vector */
130: VecGetSize(global,&n);
132: /* Set the index sets */
133: PetscMalloc(n*sizeof(PetscInt),&idx);
134: for(i=0; i<n; i++) idx[i]=i;
135:
136: /* Create local sequential vectors */
137: VecCreateSeq(PETSC_COMM_SELF,n,&tmp_vec);
139: /* Create scatter context */
140: ISCreateGeneral(PETSC_COMM_SELF,n,idx,&from);
141: ISCreateGeneral(PETSC_COMM_SELF,n,idx,&to);
142: VecScatterCreate(global,from,tmp_vec,to,&scatter);
143: VecScatterBegin(global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD,scatter);
144: VecScatterEnd(global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD,scatter);
146: VecGetArray(tmp_vec,&tmp);
147: PetscPrintf(PETSC_COMM_WORLD,"At t =%14.6e u = %14.6e %14.6e %14.6e \n",
148: time,PetscRealPart(tmp[0]),PetscRealPart(tmp[1]),PetscRealPart(tmp[2]));
149: PetscPrintf(PETSC_COMM_WORLD,"At t =%14.6e errors = %14.6e %14.6e %14.6e \n",
150: time,PetscRealPart(tmp[0]-solx(time)),PetscRealPart(tmp[1]-soly(time)),PetscRealPart(tmp[2]-solz(time)));
151: VecRestoreArray(tmp_vec,&tmp);
152: VecScatterDestroy(scatter);
153: ISDestroy(from);
154: ISDestroy(to);
155: PetscFree(idx);
156: VecDestroy(tmp_vec);
157: return 0;
158: }
162: PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec globalin,Vec globalout,void *ctx)
163: {
164: PetscScalar *inptr,*outptr;
165: PetscInt i,n,*idx;
167: IS from,to;
168: VecScatter scatter;
169: Vec tmp_in,tmp_out;
171: /* Get the length of parallel vector */
172: VecGetSize(globalin,&n);
174: /* Set the index sets */
175: PetscMalloc(n*sizeof(PetscInt),&idx);
176: for(i=0; i<n; i++) idx[i]=i;
177:
178: /* Create local sequential vectors */
179: VecCreateSeq(PETSC_COMM_SELF,n,&tmp_in);
180: VecDuplicate(tmp_in,&tmp_out);
182: /* Create scatter context */
183: ISCreateGeneral(PETSC_COMM_SELF,n,idx,&from);
184: ISCreateGeneral(PETSC_COMM_SELF,n,idx,&to);
185: VecScatterCreate(globalin,from,tmp_in,to,&scatter);
186: VecScatterBegin(globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD,scatter);
187: VecScatterEnd(globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD,scatter);
188: VecScatterDestroy(scatter);
190: /*Extract income array */
191: VecGetArray(tmp_in,&inptr);
193: /* Extract outcome array*/
194: VecGetArray(tmp_out,&outptr);
196: outptr[0] = 2.0*inptr[0]+inptr[1];
197: outptr[1] = inptr[0]+2.0*inptr[1]+inptr[2];
198: outptr[2] = inptr[1]+2.0*inptr[2];
200: VecRestoreArray(tmp_in,&inptr);
201: VecRestoreArray(tmp_out,&outptr);
203: VecScatterCreate(tmp_out,from,globalout,to,&scatter);
204: VecScatterBegin(tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD,scatter);
205: VecScatterEnd(tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD,scatter);
207: /* Destroy idx aand scatter */
208: ISDestroy(from);
209: ISDestroy(to);
210: VecScatterDestroy(scatter);
211: VecDestroy(tmp_in);
212: VecDestroy(tmp_out);
213: PetscFree(idx);
214: return 0;
215: }
219: PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec x,Mat *AA,Mat *BB,MatStructure *str,void *ctx)
220: {
221: Mat A = *AA;
222: PetscScalar v[3],*tmp;
223: PetscInt idx[3],i;
225:
226: *str = SAME_NONZERO_PATTERN;
228: idx[0]=0; idx[1]=1; idx[2]=2;
229: VecGetArray(x,&tmp);
231: i = 0;
232: v[0] = 2.0; v[1] = 1.0; v[2] = 0.0;
233: MatSetValues(A,1,&i,3,idx,v,INSERT_VALUES);
235: i = 1;
236: v[0] = 1.0; v[1] = 2.0; v[2] = 1.0;
237: MatSetValues(A,1,&i,3,idx,v,INSERT_VALUES);
238:
239: i = 2;
240: v[0]= 0.0; v[1] = 1.0; v[2] = 2.0;
241: MatSetValues(A,1,&i,3,idx,v,INSERT_VALUES);
243: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
244: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
246: VecRestoreArray(x,&tmp);
248: return 0;
249: }
251: /*
252: The exact solutions
253: */
254: PetscReal solx(PetscReal t)
255: {
256: return exp((2.0 - sqrt(2.0))*t)/2.0 - exp((2.0 - sqrt(2.0))*t)/(2.0*sqrt(2.0)) +
257: exp((2.0 + sqrt(2.0))*t)/2.0 + exp((2.0 + sqrt(2.0))*t)/(2.0*sqrt(2.0));
258: }
260: PetscReal soly(PetscReal t)
261: {
262: return exp((2.0 - sqrt(2.0))*t)/2.0 - exp((2.0 - sqrt(2.0))*t)/sqrt(2.0) +
263: exp((2.0 + sqrt(2.0))*t)/2.0 + exp((2.0 + sqrt(2.0))*t)/sqrt(2.0);
264: }
265:
266: PetscReal solz(PetscReal t)
267: {
268: return exp((2.0 - sqrt(2.0))*t)/2.0 - exp((2.0 - sqrt(2.0))*t)/(2.0*sqrt(2.0)) +
269: exp((2.0 + sqrt(2.0))*t)/2.0 + exp((2.0 + sqrt(2.0))*t)/(2.0*sqrt(2.0));
270: }