1 /* 2 Formatted test for TS routines. 3 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. 12 13 */ 14 15 static char help[] = "Solves a linear ODE. \n\n"; 16 17 #include <petscts.h> 18 #include <petscpc.h> 19 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); 25 26 extern PetscReal solx(PetscReal); 27 extern PetscReal soly(PetscReal); 28 extern PetscReal solz(PetscReal); 29 30 int main(int argc,char **argv) 31 { 32 PetscErrorCode ierr; 33 PetscInt time_steps = 100,steps; 34 PetscMPIInt size; 35 Vec global; 36 PetscReal dt,ftime; 37 TS ts; 38 Mat A = 0,S; 39 40 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 41 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr); 42 43 ierr = PetscOptionsGetInt(NULL,NULL,"-time",&time_steps,NULL);CHKERRQ(ierr); 44 45 /* set initial conditions */ 46 ierr = VecCreate(PETSC_COMM_WORLD,&global);CHKERRQ(ierr); 47 ierr = VecSetSizes(global,PETSC_DECIDE,3);CHKERRQ(ierr); 48 ierr = VecSetFromOptions(global);CHKERRQ(ierr); 49 ierr = Initial(global,NULL);CHKERRQ(ierr); 50 51 /* make timestep context */ 52 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr); 53 ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr); 54 ierr = TSMonitorSet(ts,Monitor,NULL,NULL);CHKERRQ(ierr); 55 56 dt = 0.001; 57 58 /* 59 The user provides the RHS and Jacobian 60 */ 61 ierr = TSSetRHSFunction(ts,NULL,RHSFunction,NULL);CHKERRQ(ierr); 62 ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); 63 ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,3,3);CHKERRQ(ierr); 64 ierr = MatSetFromOptions(A);CHKERRQ(ierr); 65 ierr = MatSetUp(A);CHKERRQ(ierr); 66 ierr = RHSJacobian(ts,0.0,global,A,A,NULL);CHKERRQ(ierr); 67 ierr = TSSetRHSJacobian(ts,A,A,RHSJacobian,NULL);CHKERRQ(ierr); 68 69 ierr = MatCreateShell(PETSC_COMM_WORLD,3,3,3,3,NULL,&S);CHKERRQ(ierr); 70 ierr = MatShellSetOperation(S,MATOP_MULT,(void (*)(void))MyMatMult);CHKERRQ(ierr); 71 ierr = TSSetRHSJacobian(ts,S,A,RHSJacobian,NULL);CHKERRQ(ierr); 72 73 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);CHKERRQ(ierr); 74 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 75 76 ierr = TSSetTimeStep(ts,dt);CHKERRQ(ierr); 77 ierr = TSSetMaxSteps(ts,time_steps);CHKERRQ(ierr); 78 ierr = TSSetMaxTime(ts,1);CHKERRQ(ierr); 79 ierr = TSSetSolution(ts,global);CHKERRQ(ierr); 80 81 ierr = TSSolve(ts,global);CHKERRQ(ierr); 82 ierr = TSGetSolveTime(ts,&ftime);CHKERRQ(ierr); 83 ierr = TSGetStepNumber(ts,&steps);CHKERRQ(ierr); 84 85 86 /* free the memories */ 87 88 ierr = TSDestroy(&ts);CHKERRQ(ierr); 89 ierr = VecDestroy(&global);CHKERRQ(ierr); 90 ierr = MatDestroy(&A);CHKERRQ(ierr); 91 ierr = MatDestroy(&S);CHKERRQ(ierr); 92 93 ierr = PetscFinalize(); 94 return ierr; 95 } 96 97 PetscErrorCode MyMatMult(Mat S,Vec x,Vec y) 98 { 99 PetscErrorCode ierr; 100 const PetscScalar *inptr; 101 PetscScalar *outptr; 102 103 PetscFunctionBegin; 104 ierr = VecGetArrayRead(x,&inptr);CHKERRQ(ierr); 105 ierr = VecGetArray(y,&outptr);CHKERRQ(ierr); 106 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]; 110 111 ierr = VecRestoreArrayRead(x,&inptr);CHKERRQ(ierr); 112 ierr = VecRestoreArray(y,&outptr);CHKERRQ(ierr); 113 PetscFunctionReturn(0); 114 } 115 116 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; 123 PetscErrorCode ierr; 124 125 /* determine starting point of each processor */ 126 ierr = VecGetOwnershipRange(global,&mybase,&myend);CHKERRQ(ierr); 127 ierr = VecGetLocalSize(global,&locsize);CHKERRQ(ierr); 128 129 /* Initialize the array */ 130 ierr = VecGetArray(global,&localptr);CHKERRQ(ierr); 131 for (i=0; i<locsize; i++) localptr[i] = 1.0; 132 133 if (mybase == 0) localptr[0]=1.0; 134 135 ierr = VecRestoreArray(global,&localptr);CHKERRQ(ierr); 136 return 0; 137 } 138 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; 147 148 /* Get the size of the vector */ 149 ierr = VecGetSize(global,&n);CHKERRQ(ierr); 150 151 /* Set the index sets */ 152 ierr = PetscMalloc1(n,&idx);CHKERRQ(ierr); 153 for (i=0; i<n; i++) idx[i]=i; 154 155 /* Create local sequential vectors */ 156 ierr = VecCreateSeq(PETSC_COMM_SELF,n,&tmp_vec);CHKERRQ(ierr); 157 158 /* Create scatter context */ 159 ierr = ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&from);CHKERRQ(ierr); 160 ierr = ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&to);CHKERRQ(ierr); 161 ierr = VecScatterCreate(global,from,tmp_vec,to,&scatter);CHKERRQ(ierr); 162 ierr = VecScatterBegin(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 163 ierr = VecScatterEnd(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 164 165 ierr = VecGetArrayRead(tmp_vec,&tmp);CHKERRQ(ierr); 166 ierr = 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]));CHKERRQ(ierr); 168 ierr = 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)));CHKERRQ(ierr); 170 ierr = VecRestoreArrayRead(tmp_vec,&tmp);CHKERRQ(ierr); 171 ierr = VecScatterDestroy(&scatter);CHKERRQ(ierr); 172 ierr = ISDestroy(&from);CHKERRQ(ierr); 173 ierr = ISDestroy(&to);CHKERRQ(ierr); 174 ierr = PetscFree(idx);CHKERRQ(ierr); 175 ierr = VecDestroy(&tmp_vec);CHKERRQ(ierr); 176 return 0; 177 } 178 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; 188 189 /* Get the length of parallel vector */ 190 ierr = VecGetSize(globalin,&n);CHKERRQ(ierr); 191 192 /* Set the index sets */ 193 ierr = PetscMalloc1(n,&idx);CHKERRQ(ierr); 194 for (i=0; i<n; i++) idx[i]=i; 195 196 /* Create local sequential vectors */ 197 ierr = VecCreateSeq(PETSC_COMM_SELF,n,&tmp_in);CHKERRQ(ierr); 198 ierr = VecDuplicate(tmp_in,&tmp_out);CHKERRQ(ierr); 199 200 /* Create scatter context */ 201 ierr = ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&from);CHKERRQ(ierr); 202 ierr = ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&to);CHKERRQ(ierr); 203 ierr = VecScatterCreate(globalin,from,tmp_in,to,&scatter);CHKERRQ(ierr); 204 ierr = VecScatterBegin(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 205 ierr = VecScatterEnd(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 206 ierr = VecScatterDestroy(&scatter);CHKERRQ(ierr); 207 208 /*Extract income array */ 209 ierr = VecGetArrayRead(tmp_in,&inptr);CHKERRQ(ierr); 210 211 /* Extract outcome array*/ 212 ierr = VecGetArray(tmp_out,&outptr);CHKERRQ(ierr); 213 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]; 217 218 ierr = VecRestoreArrayRead(tmp_in,&inptr);CHKERRQ(ierr); 219 ierr = VecRestoreArray(tmp_out,&outptr);CHKERRQ(ierr); 220 221 ierr = VecScatterCreate(tmp_out,from,globalout,to,&scatter);CHKERRQ(ierr); 222 ierr = VecScatterBegin(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 223 ierr = VecScatterEnd(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 224 225 /* Destroy idx aand scatter */ 226 ierr = ISDestroy(&from);CHKERRQ(ierr); 227 ierr = ISDestroy(&to);CHKERRQ(ierr); 228 ierr = VecScatterDestroy(&scatter);CHKERRQ(ierr); 229 ierr = VecDestroy(&tmp_in);CHKERRQ(ierr); 230 ierr = VecDestroy(&tmp_out);CHKERRQ(ierr); 231 ierr = PetscFree(idx);CHKERRQ(ierr); 232 return 0; 233 } 234 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; 241 242 idx[0]=0; idx[1]=1; idx[2]=2; 243 ierr = VecGetArrayRead(x,&tmp);CHKERRQ(ierr); 244 245 i = 0; 246 v[0] = 2.0; v[1] = 1.0; v[2] = 0.0; 247 ierr = MatSetValues(BB,1,&i,3,idx,v,INSERT_VALUES);CHKERRQ(ierr); 248 249 i = 1; 250 v[0] = 1.0; v[1] = 2.0; v[2] = 1.0; 251 ierr = MatSetValues(BB,1,&i,3,idx,v,INSERT_VALUES);CHKERRQ(ierr); 252 253 i = 2; 254 v[0] = 0.0; v[1] = 1.0; v[2] = 2.0; 255 ierr = MatSetValues(BB,1,&i,3,idx,v,INSERT_VALUES);CHKERRQ(ierr); 256 257 ierr = MatAssemblyBegin(BB,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 258 ierr = MatAssemblyEnd(BB,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 259 260 if (A != BB) { 261 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 262 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 263 } 264 ierr = VecRestoreArrayRead(x,&tmp);CHKERRQ(ierr); 265 266 return 0; 267 } 268 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 } 277 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 } 283 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 } 289 290 291 /*TEST 292 293 test: 294 suffix: euler 295 args: -ts_type euler 296 requires: !single 297 298 test: 299 suffix: beuler 300 args: -ts_type beuler 301 requires: !single 302 303 TEST*/ 304 305