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 petscts.h
18: #include petscpc.h
31: int main(int argc,char **argv)
32: {
34: PetscInt time_steps = 100,steps;
35: PetscMPIInt size;
36: Vec global;
37: PetscReal dt,ftime;
38: TS ts;
39: MatStructure A_structure;
40: Mat A = 0;
41:
42: PetscInitialize(&argc,&argv,(char*)0,help);
43: MPI_Comm_size(PETSC_COMM_WORLD,&size);
44:
45: PetscOptionsGetInt(PETSC_NULL,"-time",&time_steps,PETSC_NULL);
46:
47: /* set initial conditions */
48: VecCreate(PETSC_COMM_WORLD,&global);
49: VecSetSizes(global,PETSC_DECIDE,3);
50: VecSetFromOptions(global);
51: Initial(global,PETSC_NULL);
52:
53: /* make timestep context */
54: TSCreate(PETSC_COMM_WORLD,&ts);
55: TSSetProblemType(ts,TS_NONLINEAR);
56: TSMonitorSet(ts,Monitor,PETSC_NULL,PETSC_NULL);
58: dt = 0.1;
60: /*
61: The user provides the RHS and Jacobian
62: */
63: TSSetRHSFunction(ts,RHSFunction,NULL);
64: MatCreate(PETSC_COMM_WORLD,&A);
65: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,3,3);
66: MatSetFromOptions(A);
67: RHSJacobian(ts,0.0,global,&A,&A,&A_structure,NULL);
68: TSSetRHSJacobian(ts,A,A,RHSJacobian,NULL);
69:
70: TSSetFromOptions(ts);
72: TSSetInitialTimeStep(ts,0.0,dt);
73: TSSetDuration(ts,time_steps,1);
74: TSSetSolution(ts,global);
77: TSSetUp(ts);
78: TSStep(ts,&steps,&ftime);
81: /* free the memories */
82:
83: TSDestroy(ts);
84: VecDestroy(global);
85: ierr= MatDestroy(A);
87: PetscFinalize();
88: return 0;
89: }
91: /* -------------------------------------------------------------------*/
94: /* this test problem has initial values (1,1,1). */
95: PetscErrorCode Initial(Vec global,void *ctx)
96: {
97: PetscScalar *localptr;
98: PetscInt i,mybase,myend,locsize;
101: /* determine starting point of each processor */
102: VecGetOwnershipRange(global,&mybase,&myend);
103: VecGetLocalSize(global,&locsize);
105: /* Initialize the array */
106: VecGetArray(global,&localptr);
107: for (i=0; i<locsize; i++) {
108: localptr[i] = 1.0;
109: }
110:
111: if (mybase == 0) localptr[0]=1.0;
113: VecRestoreArray(global,&localptr);
114: return 0;
115: }
119: PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal time,Vec global,void *ctx)
120: {
121: VecScatter scatter;
122: IS from,to;
123: PetscInt i,n,*idx;
124: Vec tmp_vec;
126: PetscScalar *tmp;
128: /* Get the size of the vector */
129: VecGetSize(global,&n);
131: /* Set the index sets */
132: PetscMalloc(n*sizeof(PetscInt),&idx);
133: for(i=0; i<n; i++) idx[i]=i;
134:
135: /* Create local sequential vectors */
136: VecCreateSeq(PETSC_COMM_SELF,n,&tmp_vec);
138: /* Create scatter context */
139: ISCreateGeneral(PETSC_COMM_SELF,n,idx,&from);
140: ISCreateGeneral(PETSC_COMM_SELF,n,idx,&to);
141: VecScatterCreate(global,from,tmp_vec,to,&scatter);
142: VecScatterBegin(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD);
143: VecScatterEnd(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD);
145: VecGetArray(tmp_vec,&tmp);
146: PetscPrintf(PETSC_COMM_WORLD,"At t =%14.6e u = %14.6e %14.6e %14.6e \n",
147: time,PetscRealPart(tmp[0]),PetscRealPart(tmp[1]),PetscRealPart(tmp[2]));
148: PetscPrintf(PETSC_COMM_WORLD,"At t =%14.6e errors = %14.6e %14.6e %14.6e \n",
149: time,PetscRealPart(tmp[0]-solx(time)),PetscRealPart(tmp[1]-soly(time)),PetscRealPart(tmp[2]-solz(time)));
150: VecRestoreArray(tmp_vec,&tmp);
151: VecScatterDestroy(scatter);
152: ISDestroy(from);
153: ISDestroy(to);
154: PetscFree(idx);
155: VecDestroy(tmp_vec);
156: return 0;
157: }
161: PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec globalin,Vec globalout,void *ctx)
162: {
163: PetscScalar *inptr,*outptr;
164: PetscInt i,n,*idx;
166: IS from,to;
167: VecScatter scatter;
168: Vec tmp_in,tmp_out;
170: /* Get the length of parallel vector */
171: VecGetSize(globalin,&n);
173: /* Set the index sets */
174: PetscMalloc(n*sizeof(PetscInt),&idx);
175: for(i=0; i<n; i++) idx[i]=i;
176:
177: /* Create local sequential vectors */
178: VecCreateSeq(PETSC_COMM_SELF,n,&tmp_in);
179: VecDuplicate(tmp_in,&tmp_out);
181: /* Create scatter context */
182: ISCreateGeneral(PETSC_COMM_SELF,n,idx,&from);
183: ISCreateGeneral(PETSC_COMM_SELF,n,idx,&to);
184: VecScatterCreate(globalin,from,tmp_in,to,&scatter);
185: VecScatterBegin(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD);
186: VecScatterEnd(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD);
187: VecScatterDestroy(scatter);
189: /*Extract income array */
190: VecGetArray(tmp_in,&inptr);
192: /* Extract outcome array*/
193: VecGetArray(tmp_out,&outptr);
195: outptr[0] = 2.0*inptr[0]+inptr[1];
196: outptr[1] = inptr[0]+2.0*inptr[1]+inptr[2];
197: outptr[2] = inptr[1]+2.0*inptr[2];
199: VecRestoreArray(tmp_in,&inptr);
200: VecRestoreArray(tmp_out,&outptr);
202: VecScatterCreate(tmp_out,from,globalout,to,&scatter);
203: VecScatterBegin(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD);
204: VecScatterEnd(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD);
206: /* Destroy idx aand scatter */
207: ISDestroy(from);
208: ISDestroy(to);
209: VecScatterDestroy(scatter);
210: VecDestroy(tmp_in);
211: VecDestroy(tmp_out);
212: PetscFree(idx);
213: return 0;
214: }
218: PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec x,Mat *AA,Mat *BB,MatStructure *str,void *ctx)
219: {
220: Mat A = *AA;
221: PetscScalar v[3],*tmp;
222: PetscInt idx[3],i;
224:
225: *str = SAME_NONZERO_PATTERN;
227: idx[0]=0; idx[1]=1; idx[2]=2;
228: VecGetArray(x,&tmp);
230: i = 0;
231: v[0] = 2.0; v[1] = 1.0; v[2] = 0.0;
232: MatSetValues(A,1,&i,3,idx,v,INSERT_VALUES);
234: i = 1;
235: v[0] = 1.0; v[1] = 2.0; v[2] = 1.0;
236: MatSetValues(A,1,&i,3,idx,v,INSERT_VALUES);
237:
238: i = 2;
239: v[0]= 0.0; v[1] = 1.0; v[2] = 2.0;
240: MatSetValues(A,1,&i,3,idx,v,INSERT_VALUES);
242: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
243: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
245: VecRestoreArray(x,&tmp);
247: return 0;
248: }
250: /*
251: The exact solutions
252: */
253: PetscReal solx(PetscReal t)
254: {
255: return exp((2.0 - sqrt(2.0))*t)/2.0 - exp((2.0 - sqrt(2.0))*t)/(2.0*sqrt(2.0)) +
256: exp((2.0 + sqrt(2.0))*t)/2.0 + exp((2.0 + sqrt(2.0))*t)/(2.0*sqrt(2.0));
257: }
259: PetscReal soly(PetscReal t)
260: {
261: return exp((2.0 - sqrt(2.0))*t)/2.0 - exp((2.0 - sqrt(2.0))*t)/sqrt(2.0) +
262: exp((2.0 + sqrt(2.0))*t)/2.0 + exp((2.0 + sqrt(2.0))*t)/sqrt(2.0);
263: }
264:
265: PetscReal solz(PetscReal t)
266: {
267: return exp((2.0 - sqrt(2.0))*t)/2.0 - exp((2.0 - sqrt(2.0))*t)/(2.0*sqrt(2.0)) +
268: exp((2.0 + sqrt(2.0))*t)/2.0 + exp((2.0 + sqrt(2.0))*t)/(2.0*sqrt(2.0));
269: }