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 Vec global; 35 PetscReal dt,ftime; 36 TS ts; 37 Mat A = 0,S; 38 39 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 40 ierr = PetscOptionsGetInt(NULL,NULL,"-time",&time_steps,NULL);CHKERRQ(ierr); 41 42 /* set initial conditions */ 43 ierr = VecCreate(PETSC_COMM_WORLD,&global);CHKERRQ(ierr); 44 ierr = VecSetSizes(global,PETSC_DECIDE,3);CHKERRQ(ierr); 45 ierr = VecSetFromOptions(global);CHKERRQ(ierr); 46 ierr = Initial(global,NULL);CHKERRQ(ierr); 47 48 /* make timestep context */ 49 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr); 50 ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr); 51 ierr = TSMonitorSet(ts,Monitor,NULL,NULL);CHKERRQ(ierr); 52 dt = 0.001; 53 54 /* 55 The user provides the RHS and Jacobian 56 */ 57 ierr = TSSetRHSFunction(ts,NULL,RHSFunction,NULL);CHKERRQ(ierr); 58 ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); 59 ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,3,3);CHKERRQ(ierr); 60 ierr = MatSetFromOptions(A);CHKERRQ(ierr); 61 ierr = MatSetUp(A);CHKERRQ(ierr); 62 ierr = RHSJacobian(ts,0.0,global,A,A,NULL);CHKERRQ(ierr); 63 ierr = TSSetRHSJacobian(ts,A,A,RHSJacobian,NULL);CHKERRQ(ierr); 64 65 ierr = MatCreateShell(PETSC_COMM_WORLD,3,3,3,3,NULL,&S);CHKERRQ(ierr); 66 ierr = MatShellSetOperation(S,MATOP_MULT,(void (*)(void))MyMatMult);CHKERRQ(ierr); 67 ierr = TSSetRHSJacobian(ts,S,A,RHSJacobian,NULL);CHKERRQ(ierr); 68 69 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);CHKERRQ(ierr); 70 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 71 72 ierr = TSSetTimeStep(ts,dt);CHKERRQ(ierr); 73 ierr = TSSetMaxSteps(ts,time_steps);CHKERRQ(ierr); 74 ierr = TSSetMaxTime(ts,1);CHKERRQ(ierr); 75 ierr = TSSetSolution(ts,global);CHKERRQ(ierr); 76 77 ierr = TSSolve(ts,global);CHKERRQ(ierr); 78 ierr = TSGetSolveTime(ts,&ftime);CHKERRQ(ierr); 79 ierr = TSGetStepNumber(ts,&steps);CHKERRQ(ierr); 80 81 /* free the memories */ 82 83 ierr = TSDestroy(&ts);CHKERRQ(ierr); 84 ierr = VecDestroy(&global);CHKERRQ(ierr); 85 ierr = MatDestroy(&A);CHKERRQ(ierr); 86 ierr = MatDestroy(&S);CHKERRQ(ierr); 87 88 ierr = PetscFinalize(); 89 return ierr; 90 } 91 92 PetscErrorCode MyMatMult(Mat S,Vec x,Vec y) 93 { 94 PetscErrorCode ierr; 95 const PetscScalar *inptr; 96 PetscScalar *outptr; 97 98 PetscFunctionBegin; 99 ierr = VecGetArrayRead(x,&inptr);CHKERRQ(ierr); 100 ierr = VecGetArrayWrite(y,&outptr);CHKERRQ(ierr); 101 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]; 105 106 ierr = VecRestoreArrayRead(x,&inptr);CHKERRQ(ierr); 107 ierr = VecRestoreArrayWrite(y,&outptr);CHKERRQ(ierr); 108 PetscFunctionReturn(0); 109 } 110 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; 116 PetscErrorCode ierr; 117 118 /* determine starting point of each processor */ 119 ierr = VecGetOwnershipRange(global,&mybase,&myend);CHKERRQ(ierr); 120 ierr = VecGetLocalSize(global,&locsize);CHKERRQ(ierr); 121 122 /* Initialize the array */ 123 ierr = VecGetArrayWrite(global,&localptr);CHKERRQ(ierr); 124 for (i=0; i<locsize; i++) localptr[i] = 1.0; 125 126 if (mybase == 0) localptr[0]=1.0; 127 128 ierr = VecRestoreArrayWrite(global,&localptr);CHKERRQ(ierr); 129 return 0; 130 } 131 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; 140 141 /* Get the size of the vector */ 142 ierr = VecGetSize(global,&n);CHKERRQ(ierr); 143 144 /* Set the index sets */ 145 ierr = PetscMalloc1(n,&idx);CHKERRQ(ierr); 146 for (i=0; i<n; i++) idx[i]=i; 147 148 /* Create local sequential vectors */ 149 ierr = VecCreateSeq(PETSC_COMM_SELF,n,&tmp_vec);CHKERRQ(ierr); 150 151 /* Create scatter context */ 152 ierr = ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&from);CHKERRQ(ierr); 153 ierr = ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&to);CHKERRQ(ierr); 154 ierr = VecScatterCreate(global,from,tmp_vec,to,&scatter);CHKERRQ(ierr); 155 ierr = VecScatterBegin(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 156 ierr = VecScatterEnd(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 157 158 ierr = VecGetArrayRead(tmp_vec,&tmp);CHKERRQ(ierr); 159 ierr = 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]));CHKERRQ(ierr); 161 ierr = 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)));CHKERRQ(ierr); 163 ierr = VecRestoreArrayRead(tmp_vec,&tmp);CHKERRQ(ierr); 164 ierr = VecScatterDestroy(&scatter);CHKERRQ(ierr); 165 ierr = ISDestroy(&from);CHKERRQ(ierr); 166 ierr = ISDestroy(&to);CHKERRQ(ierr); 167 ierr = PetscFree(idx);CHKERRQ(ierr); 168 ierr = VecDestroy(&tmp_vec);CHKERRQ(ierr); 169 return 0; 170 } 171 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; 181 182 /* Get the length of parallel vector */ 183 ierr = VecGetSize(globalin,&n);CHKERRQ(ierr); 184 185 /* Set the index sets */ 186 ierr = PetscMalloc1(n,&idx);CHKERRQ(ierr); 187 for (i=0; i<n; i++) idx[i]=i; 188 189 /* Create local sequential vectors */ 190 ierr = VecCreateSeq(PETSC_COMM_SELF,n,&tmp_in);CHKERRQ(ierr); 191 ierr = VecDuplicate(tmp_in,&tmp_out);CHKERRQ(ierr); 192 193 /* Create scatter context */ 194 ierr = ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&from);CHKERRQ(ierr); 195 ierr = ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&to);CHKERRQ(ierr); 196 ierr = VecScatterCreate(globalin,from,tmp_in,to,&scatter);CHKERRQ(ierr); 197 ierr = VecScatterBegin(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 198 ierr = VecScatterEnd(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 199 ierr = VecScatterDestroy(&scatter);CHKERRQ(ierr); 200 201 /*Extract income array */ 202 ierr = VecGetArrayRead(tmp_in,&inptr);CHKERRQ(ierr); 203 204 /* Extract outcome array*/ 205 ierr = VecGetArrayWrite(tmp_out,&outptr);CHKERRQ(ierr); 206 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]; 210 211 ierr = VecRestoreArrayRead(tmp_in,&inptr);CHKERRQ(ierr); 212 ierr = VecRestoreArrayWrite(tmp_out,&outptr);CHKERRQ(ierr); 213 214 ierr = VecScatterCreate(tmp_out,from,globalout,to,&scatter);CHKERRQ(ierr); 215 ierr = VecScatterBegin(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 216 ierr = VecScatterEnd(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 217 218 /* Destroy idx aand scatter */ 219 ierr = ISDestroy(&from);CHKERRQ(ierr); 220 ierr = ISDestroy(&to);CHKERRQ(ierr); 221 ierr = VecScatterDestroy(&scatter);CHKERRQ(ierr); 222 ierr = VecDestroy(&tmp_in);CHKERRQ(ierr); 223 ierr = VecDestroy(&tmp_out);CHKERRQ(ierr); 224 ierr = PetscFree(idx);CHKERRQ(ierr); 225 return 0; 226 } 227 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; 234 235 idx[0]=0; idx[1]=1; idx[2]=2; 236 ierr = VecGetArrayRead(x,&tmp);CHKERRQ(ierr); 237 238 i = 0; 239 v[0] = 2.0; v[1] = 1.0; v[2] = 0.0; 240 ierr = MatSetValues(BB,1,&i,3,idx,v,INSERT_VALUES);CHKERRQ(ierr); 241 242 i = 1; 243 v[0] = 1.0; v[1] = 2.0; v[2] = 1.0; 244 ierr = MatSetValues(BB,1,&i,3,idx,v,INSERT_VALUES);CHKERRQ(ierr); 245 246 i = 2; 247 v[0] = 0.0; v[1] = 1.0; v[2] = 2.0; 248 ierr = MatSetValues(BB,1,&i,3,idx,v,INSERT_VALUES);CHKERRQ(ierr); 249 250 ierr = MatAssemblyBegin(BB,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 251 ierr = MatAssemblyEnd(BB,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 252 253 if (A != BB) { 254 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 255 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 256 } 257 ierr = VecRestoreArrayRead(x,&tmp);CHKERRQ(ierr); 258 259 return 0; 260 } 261 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 } 270 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 } 276 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 } 282 283 /*TEST 284 285 test: 286 suffix: euler 287 args: -ts_type euler 288 requires: !single 289 290 test: 291 suffix: beuler 292 args: -ts_type beuler 293 requires: !single 294 295 TEST*/ 296 297