Actual source code: ex2.c
1: /*$Id: ex2.c,v 1.36 2001/08/07 21:31:27 bsmith Exp $*/
2: /*
3: Formatted test for TS routines.
5: Solves U_t=F(t,u)
6: Where:
7:
8: [2*u1+u2
9: F(t,u)= [u1+2*u2+u3
10: [ u2+2*u3
11: We can compare the solutions from euler, beuler and PVODE to
12: see what is the difference.
14: */
16: static char help[] = "Solves a nonlinear ODE. \n\n";
18: #include petscsys.h
19: #include petscts.h
20: #include petscpc.h
22: extern int RHSFunction(TS,PetscReal,Vec,Vec,void*);
23: extern int RHSJacobian(TS,PetscReal,Vec,Mat*,Mat*,MatStructure *,void*);
24: extern int Monitor(TS,int,PetscReal,Vec,void *);
25: extern int Initial(Vec,void *);
27: extern PetscReal solx(PetscReal);
28: extern PetscReal soly(PetscReal);
29: extern PetscReal solz(PetscReal);
33: int main(int argc,char **argv)
34: {
35: int ierr,time_steps = 100,steps,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: TSSetMonitor(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,PETSC_DECIDE,PETSC_DECIDE,3,3,&A);
65: MatSetFromOptions(A);
66: RHSJacobian(ts,0.0,global,&A,&A,&A_structure,NULL);
67: TSSetRHSJacobian(ts,A,A,RHSJacobian,NULL);
68:
69: TSSetFromOptions(ts);
71: TSSetInitialTimeStep(ts,0.0,dt);
72: TSSetDuration(ts,time_steps,1);
73: TSSetSolution(ts,global);
76: TSSetUp(ts);
77: TSStep(ts,&steps,&ftime);
80: /* free the memories */
81:
82: TSDestroy(ts);
83: VecDestroy(global);
84: ierr= MatDestroy(A);
86: PetscFinalize();
87: return 0;
88: }
90: /* -------------------------------------------------------------------*/
93: /* this test problem has initial values (1,1,1). */
94: int Initial(Vec global,void *ctx)
95: {
96: PetscScalar *localptr;
97: int i,mybase,myend,ierr,locsize;
99: /* determine starting point of each processor */
100: VecGetOwnershipRange(global,&mybase,&myend);
101: VecGetLocalSize(global,&locsize);
103: /* Initialize the array */
104: VecGetArray(global,&localptr);
105: for (i=0; i<locsize; i++) {
106: localptr[i] = 1.0;
107: }
108:
109: if (mybase == 0) localptr[0]=1.0;
111: VecRestoreArray(global,&localptr);
112: return 0;
113: }
117: int Monitor(TS ts,int step,PetscReal time,Vec global,void *ctx)
118: {
119: VecScatter scatter;
120: IS from,to;
121: int i,n,*idx;
122: Vec tmp_vec;
123: int ierr;
124: PetscScalar *tmp;
126: /* Get the size of the vector */
127: VecGetSize(global,&n);
129: /* Set the index sets */
130: PetscMalloc(n*sizeof(int),&idx);
131: for(i=0; i<n; i++) idx[i]=i;
132:
133: /* Create local sequential vectors */
134: VecCreateSeq(PETSC_COMM_SELF,n,&tmp_vec);
136: /* Create scatter context */
137: ISCreateGeneral(PETSC_COMM_SELF,n,idx,&from);
138: ISCreateGeneral(PETSC_COMM_SELF,n,idx,&to);
139: VecScatterCreate(global,from,tmp_vec,to,&scatter);
140: VecScatterBegin(global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD,scatter);
141: VecScatterEnd(global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD,scatter);
143: VecGetArray(tmp_vec,&tmp);
144: PetscPrintf(PETSC_COMM_WORLD,"At t =%14.6e u = %14.6e %14.6e %14.6e \n",
145: time,PetscRealPart(tmp[0]),PetscRealPart(tmp[1]),PetscRealPart(tmp[2]));
146: PetscPrintf(PETSC_COMM_WORLD,"At t =%14.6e errors = %14.6e %14.6e %14.6e \n",
147: time,PetscRealPart(tmp[0]-solx(time)),PetscRealPart(tmp[1]-soly(time)),PetscRealPart(tmp[2]-solz(time)));
148: VecRestoreArray(tmp_vec,&tmp);
149: VecScatterDestroy(scatter);
150: ISDestroy(from);
151: ISDestroy(to);
152: PetscFree(idx);
153: VecDestroy(tmp_vec);
154: return 0;
155: }
159: int RHSFunction(TS ts,PetscReal t,Vec globalin,Vec globalout,void *ctx)
160: {
161: PetscScalar *inptr,*outptr;
162: int i,n,ierr,*idx;
163: IS from,to;
164: VecScatter scatter;
165: Vec tmp_in,tmp_out;
167: /* Get the length of parallel vector */
168: VecGetSize(globalin,&n);
170: /* Set the index sets */
171: PetscMalloc(n*sizeof(int),&idx);
172: for(i=0; i<n; i++) idx[i]=i;
173:
174: /* Create local sequential vectors */
175: VecCreateSeq(PETSC_COMM_SELF,n,&tmp_in);
176: VecDuplicate(tmp_in,&tmp_out);
178: /* Create scatter context */
179: ISCreateGeneral(PETSC_COMM_SELF,n,idx,&from);
180: ISCreateGeneral(PETSC_COMM_SELF,n,idx,&to);
181: VecScatterCreate(globalin,from,tmp_in,to,&scatter);
182: VecScatterBegin(globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD,scatter);
183: VecScatterEnd(globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD,scatter);
184: VecScatterDestroy(scatter);
186: /*Extract income array */
187: VecGetArray(tmp_in,&inptr);
189: /* Extract outcome array*/
190: VecGetArray(tmp_out,&outptr);
192: outptr[0] = 2.0*inptr[0]+inptr[1];
193: outptr[1] = inptr[0]+2.0*inptr[1]+inptr[2];
194: outptr[2] = inptr[1]+2.0*inptr[2];
196: VecRestoreArray(tmp_in,&inptr);
197: VecRestoreArray(tmp_out,&outptr);
199: VecScatterCreate(tmp_out,from,globalout,to,&scatter);
200: VecScatterBegin(tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD,scatter);
201: VecScatterEnd(tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD,scatter);
203: /* Destroy idx aand scatter */
204: ISDestroy(from);
205: ISDestroy(to);
206: VecScatterDestroy(scatter);
207: VecDestroy(tmp_in);
208: VecDestroy(tmp_out);
209: PetscFree(idx);
210: return 0;
211: }
215: int RHSJacobian(TS ts,PetscReal t,Vec x,Mat *AA,Mat *BB,MatStructure *str,void *ctx)
216: {
217: Mat A = *AA;
218: PetscScalar v[3],*tmp;
219: int idx[3],i,ierr;
220:
221: *str = SAME_NONZERO_PATTERN;
223: idx[0]=0; idx[1]=1; idx[2]=2;
224: VecGetArray(x,&tmp);
226: i = 0;
227: v[0] = 2.0; v[1] = 1.0; v[2] = 0.0;
228: MatSetValues(A,1,&i,3,idx,v,INSERT_VALUES);
230: i = 1;
231: v[0] = 1.0; v[1] = 2.0; v[2] = 1.0;
232: MatSetValues(A,1,&i,3,idx,v,INSERT_VALUES);
233:
234: i = 2;
235: v[0]= 0.0; v[1] = 1.0; v[2] = 2.0;
236: MatSetValues(A,1,&i,3,idx,v,INSERT_VALUES);
238: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
239: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
241: VecRestoreArray(x,&tmp);
243: return 0;
244: }
246: /*
247: The exact solutions
248: */
249: PetscReal solx(PetscReal t)
250: {
251: return exp((2.0 - sqrt(2.0))*t)/2.0 - exp((2.0 - sqrt(2.0))*t)/(2.0*sqrt(2.0)) +
252: exp((2.0 + sqrt(2.0))*t)/2.0 + exp((2.0 + sqrt(2.0))*t)/(2.0*sqrt(2.0));
253: }
255: PetscReal soly(PetscReal t)
256: {
257: return exp((2.0 - sqrt(2.0))*t)/2.0 - exp((2.0 - sqrt(2.0))*t)/sqrt(2.0) +
258: exp((2.0 + sqrt(2.0))*t)/2.0 + exp((2.0 + sqrt(2.0))*t)/sqrt(2.0);
259: }
260:
261: PetscReal solz(PetscReal t)
262: {
263: return exp((2.0 - sqrt(2.0))*t)/2.0 - exp((2.0 - sqrt(2.0))*t)/(2.0*sqrt(2.0)) +
264: exp((2.0 + sqrt(2.0))*t)/2.0 + exp((2.0 + sqrt(2.0))*t)/(2.0*sqrt(2.0));
265: }