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