1 /* 2 Note: 3 To performance a convergence study, a reference solution can be generated with a small stepsize and stored into a binary file, e.g. -ts_type ssp -ts_dt 1e-5 -ts_max_steps 30000 -ts_view_solution binary:reference.bin 4 Errors can be computed in the following runs with -simulation -f reference.bin 5 6 Multirate RK options: 7 -ts_rk_dtratio is the ratio between larger time step size and small time step size 8 -ts_rk_multirate_type has three choices: none (for single step RK) 9 combined (for multirate method and user just need to provide the same RHS with the single step RK) 10 partitioned (for multiraet method and user need to provide two RHS, one is for fast components and the orther is for slow components 11 */ 12 13 static const char help[] = "1D periodic Finite Volume solver in slope-limiter form with semidiscrete time stepping.\n" 14 " advection - Variable coefficient scalar advection\n" 15 " u_t + (a*u)_x = 0\n" 16 " a is a piecewise function with two values say a0 and a1 (both a0 and a1 are constant).\n" 17 " in this toy problem we assume a=a0 when x<0 and a=a1 when x>0 with a0<a1 which has the same discontinuous point with initial condition.\n" 18 " we don't provide an exact solution, so you need to generate reference solution to do convergnece staudy,\n" 19 " more precisely, you need firstly generate a reference to a binary file say file.bin, then on commend line,\n" 20 " you should type -simulation -f file.bin.\n" 21 " you can choose the number of grids by -da_grid_x.\n" 22 " you can choose the value of a by -physics_advect_a1 and -physics_advect_a2.\n"; 23 24 #include <petscts.h> 25 #include <petscdm.h> 26 #include <petscdmda.h> 27 #include <petscdraw.h> 28 #include <petsc/private/tsimpl.h> 29 30 #include <petsc/private/kernels/blockinvert.h> /* For the Kernel_*_gets_* stuff for BAIJ */ 31 32 #include "finitevolume1d.h" 33 34 static inline PetscReal RangeMod(PetscReal a,PetscReal xmin,PetscReal xmax) { PetscReal range = xmax-xmin; return xmin +PetscFmodReal(range+PetscFmodReal(a,range),range); } 35 36 /* --------------------------------- Advection ----------------------------------- */ 37 38 typedef struct { 39 PetscReal a[2]; /* advective velocity */ 40 } AdvectCtx; 41 42 static PetscErrorCode PhysicsRiemann_Advect(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed,PetscReal x,PetscReal xmin, PetscReal xmax) 43 { 44 AdvectCtx *ctx = (AdvectCtx*)vctx; 45 PetscReal *speed; 46 47 PetscFunctionBeginUser; 48 speed = ctx->a; 49 if (x==0 || x == xmin || x == xmax) flux[0] = PetscMax(0,(speed[0]+speed[1])/2.0)*uL[0] + PetscMin(0,(speed[0]+speed[1])/2.0)*uR[0]; /* if condition need to be changed base on different problem, '0' is the discontinuous point of a */ 50 else if (x<0) flux[0] = PetscMax(0,speed[0])*uL[0] + PetscMin(0,speed[0])*uR[0]; /* else if condition need to be changed based on diferent problem, 'x = 0' is discontinuous point of a */ 51 else flux[0] = PetscMax(0,speed[1])*uL[0] + PetscMin(0,speed[1])*uR[0]; 52 *maxspeed = *speed; 53 PetscFunctionReturn(0); 54 } 55 56 static PetscErrorCode PhysicsCharacteristic_Advect(void *vctx,PetscInt m,const PetscScalar *u,PetscScalar *X,PetscScalar *Xi,PetscReal *speeds,PetscReal x) 57 { 58 AdvectCtx *ctx = (AdvectCtx*)vctx; 59 60 PetscFunctionBeginUser; 61 X[0] = 1.; 62 Xi[0] = 1.; 63 if (x<0) speeds[0] = ctx->a[0]; /* x is discontinuous point of a */ 64 else speeds[0] = ctx->a[1]; 65 PetscFunctionReturn(0); 66 } 67 68 static PetscErrorCode PhysicsSample_Advect(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u) 69 { 70 AdvectCtx *ctx = (AdvectCtx*)vctx; 71 PetscReal *a = ctx->a,x0; 72 73 PetscFunctionBeginUser; 74 if (x<0){ /* x is cell center */ 75 switch (bctype) { 76 case FVBC_OUTFLOW: x0 = x-a[0]*t; break; 77 case FVBC_PERIODIC: x0 = RangeMod(x-a[0]*t,xmin,xmax); break; 78 default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown BCType"); 79 } 80 switch (initial) { 81 case 0: u[0] = (x0 < 0) ? 1 : -1; break; 82 case 1: u[0] = (x0 < 0) ? -1 : 1; break; 83 case 2: u[0] = (0 < x0 && x0 < 1) ? 1 : 0; break; 84 case 3: u[0] = PetscSinReal(2*PETSC_PI*x0); break; 85 case 4: u[0] = PetscAbs(x0); break; 86 case 5: u[0] = (x0 < 0 || x0 > 0.5) ? 0 : PetscSqr(PetscSinReal(2*PETSC_PI*x0)); break; 87 case 6: u[0] = (x0 < 0) ? 0 : ((x0 < 1) ? x0 : ((x0 < 2) ? 2-x0 : 0)); break; 88 case 7: u[0] = PetscPowReal(PetscSinReal(PETSC_PI*x0),10.0);break; 89 default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown initial condition"); 90 } 91 } 92 else{ 93 switch (bctype) { 94 case FVBC_OUTFLOW: x0 = x-a[1]*t; break; 95 case FVBC_PERIODIC: x0 = RangeMod(x-a[1]*t,xmin,xmax); break; 96 default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown BCType"); 97 } 98 switch (initial) { 99 case 0: u[0] = (x0 < 0) ? 1 : -1; break; 100 case 1: u[0] = (x0 < 0) ? -1 : 1; break; 101 case 2: u[0] = (0 < x0 && x0 < 1) ? 1 : 0; break; 102 case 3: u[0] = PetscSinReal(2*PETSC_PI*x0); break; 103 case 4: u[0] = PetscAbs(x0); break; 104 case 5: u[0] = (x0 < 0 || x0 > 0.5) ? 0 : PetscSqr(PetscSinReal(2*PETSC_PI*x0)); break; 105 case 6: u[0] = (x0 < 0) ? 0 : ((x0 < 1) ? x0 : ((x0 < 2) ? 2-x0 : 0)); break; 106 case 7: u[0] = PetscPowReal(PetscSinReal(PETSC_PI*x0),10.0);break; 107 default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown initial condition"); 108 } 109 } 110 PetscFunctionReturn(0); 111 } 112 113 static PetscErrorCode PhysicsCreate_Advect(FVCtx *ctx) 114 { 115 PetscErrorCode ierr; 116 AdvectCtx *user; 117 118 PetscFunctionBeginUser; 119 CHKERRQ(PetscNew(&user)); 120 ctx->physics.sample = PhysicsSample_Advect; 121 ctx->physics.riemann = PhysicsRiemann_Advect; 122 ctx->physics.characteristic = PhysicsCharacteristic_Advect; 123 ctx->physics.destroy = PhysicsDestroy_SimpleFree; 124 ctx->physics.user = user; 125 ctx->physics.dof = 1; 126 CHKERRQ(PetscStrallocpy("u",&ctx->physics.fieldname[0])); 127 user->a[0] = 1; 128 user->a[1] = 1; 129 ierr = PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for advection","");CHKERRQ(ierr); 130 { 131 CHKERRQ(PetscOptionsReal("-physics_advect_a1","Speed1","",user->a[0],&user->a[0],NULL)); 132 CHKERRQ(PetscOptionsReal("-physics_advect_a2","Speed2","",user->a[1],&user->a[1],NULL)); 133 } 134 ierr = PetscOptionsEnd();CHKERRQ(ierr); 135 PetscFunctionReturn(0); 136 } 137 138 /* --------------------------------- Finite Volume Solver for slow parts ----------------------------------- */ 139 140 PetscErrorCode FVRHSFunctionslow(TS ts,PetscReal time,Vec X,Vec F,void *vctx) 141 { 142 FVCtx *ctx = (FVCtx*)vctx; 143 PetscInt i,j,k,Mx,dof,xs,xm,len_slow; 144 PetscReal hx,cfl_idt = 0; 145 PetscScalar *x,*f,*slope; 146 Vec Xloc; 147 DM da; 148 149 PetscFunctionBeginUser; 150 CHKERRQ(TSGetDM(ts,&da)); 151 CHKERRQ(DMGetLocalVector(da,&Xloc)); 152 CHKERRQ(DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0)); 153 hx = (ctx->xmax-ctx->xmin)/Mx; 154 CHKERRQ(DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc)); 155 CHKERRQ(DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc)); 156 CHKERRQ(VecZeroEntries(F)); 157 CHKERRQ(DMDAVecGetArray(da,Xloc,&x)); 158 CHKERRQ(VecGetArray(F,&f)); 159 CHKERRQ(DMDAGetArray(da,PETSC_TRUE,&slope)); 160 CHKERRQ(DMDAGetCorners(da,&xs,0,0,&xm,0,0)); 161 CHKERRQ(ISGetSize(ctx->iss,&len_slow)); 162 163 if (ctx->bctype == FVBC_OUTFLOW) { 164 for (i=xs-2; i<0; i++) { 165 for (j=0; j<dof; j++) x[i*dof+j] = x[j]; 166 } 167 for (i=Mx; i<xs+xm+2; i++) { 168 for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j]; 169 } 170 } 171 for (i=xs-1; i<xs+xm+1; i++) { 172 struct _LimitInfo info; 173 PetscScalar *cjmpL,*cjmpR; 174 if (i < len_slow+1) { 175 /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */ 176 CHKERRQ((*ctx->physics.characteristic)(ctx->physics.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds,ctx->xmin+hx*i)); 177 /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */ 178 CHKERRQ(PetscArrayzero(ctx->cjmpLR,2*dof)); 179 cjmpL = &ctx->cjmpLR[0]; 180 cjmpR = &ctx->cjmpLR[dof]; 181 for (j=0; j<dof; j++) { 182 PetscScalar jmpL,jmpR; 183 jmpL = x[(i+0)*dof+j]-x[(i-1)*dof+j]; 184 jmpR = x[(i+1)*dof+j]-x[(i+0)*dof+j]; 185 for (k=0; k<dof; k++) { 186 cjmpL[k] += ctx->Rinv[k+j*dof]*jmpL; 187 cjmpR[k] += ctx->Rinv[k+j*dof]*jmpR; 188 } 189 } 190 /* Apply limiter to the left and right characteristic jumps */ 191 info.m = dof; 192 info.hx = hx; 193 (*ctx->limit)(&info,cjmpL,cjmpR,ctx->cslope); 194 for (j=0; j<dof; j++) ctx->cslope[j] /= hx; /* rescale to a slope */ 195 for (j=0; j<dof; j++) { 196 PetscScalar tmp = 0; 197 for (k=0; k<dof; k++) tmp += ctx->R[j+k*dof]*ctx->cslope[k]; 198 slope[i*dof+j] = tmp; 199 } 200 } 201 } 202 203 for (i=xs; i<xs+xm+1; i++) { 204 PetscReal maxspeed; 205 PetscScalar *uL,*uR; 206 if (i < len_slow) { /* slow parts can be changed based on a */ 207 uL = &ctx->uLR[0]; 208 uR = &ctx->uLR[dof]; 209 for (j=0; j<dof; j++) { 210 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hx/2; 211 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hx/2; 212 } 213 CHKERRQ((*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed,ctx->xmin+hx*i,ctx->xmin,ctx->xmax)); 214 cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hx)); /* Max allowable value of 1/Delta t */ 215 if (i > xs) { 216 for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hx; 217 } 218 if (i < xs+xm) { 219 for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hx; 220 } 221 } 222 if (i == len_slow) { /* slow parts can be changed based on a */ 223 uL = &ctx->uLR[0]; 224 uR = &ctx->uLR[dof]; 225 for (j=0; j<dof; j++) { 226 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hx/2; 227 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hx/2; 228 } 229 CHKERRQ((*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed,ctx->xmin+hx*i,ctx->xmin,ctx->xmax)); 230 cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hx)); /* Max allowable value of 1/Delta t */ 231 if (i > xs) { 232 for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hx; 233 } 234 } 235 } 236 CHKERRQ(DMDAVecRestoreArray(da,Xloc,&x)); 237 CHKERRQ(VecRestoreArray(F,&f)); 238 CHKERRQ(DMDARestoreArray(da,PETSC_TRUE,&slope)); 239 CHKERRQ(DMRestoreLocalVector(da,&Xloc)); 240 PetscFunctionReturn(0); 241 } 242 243 /* --------------------------------- Finite Volume Solver for fast parts ----------------------------------- */ 244 245 PetscErrorCode FVRHSFunctionfast(TS ts,PetscReal time,Vec X,Vec F,void *vctx) 246 { 247 FVCtx *ctx = (FVCtx*)vctx; 248 PetscInt i,j,k,Mx,dof,xs,xm,len_slow; 249 PetscReal hx,cfl_idt = 0; 250 PetscScalar *x,*f,*slope; 251 Vec Xloc; 252 DM da; 253 254 PetscFunctionBeginUser; 255 CHKERRQ(TSGetDM(ts,&da)); 256 CHKERRQ(DMGetLocalVector(da,&Xloc)); 257 CHKERRQ(DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0)); 258 hx = (ctx->xmax-ctx->xmin)/Mx; 259 CHKERRQ(DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc)); 260 CHKERRQ(DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc)); 261 CHKERRQ(VecZeroEntries(F)); 262 CHKERRQ(DMDAVecGetArray(da,Xloc,&x)); 263 CHKERRQ(VecGetArray(F,&f)); 264 CHKERRQ(DMDAGetArray(da,PETSC_TRUE,&slope)); 265 CHKERRQ(DMDAGetCorners(da,&xs,0,0,&xm,0,0)); 266 CHKERRQ(ISGetSize(ctx->iss,&len_slow)); 267 268 if (ctx->bctype == FVBC_OUTFLOW) { 269 for (i=xs-2; i<0; i++) { 270 for (j=0; j<dof; j++) x[i*dof+j] = x[j]; 271 } 272 for (i=Mx; i<xs+xm+2; i++) { 273 for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j]; 274 } 275 } 276 for (i=xs-1; i<xs+xm+1; i++) { 277 struct _LimitInfo info; 278 PetscScalar *cjmpL,*cjmpR; 279 if (i > len_slow-2) { 280 CHKERRQ((*ctx->physics.characteristic)(ctx->physics.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds,ctx->xmin+hx*i)); 281 CHKERRQ(PetscArrayzero(ctx->cjmpLR,2*dof)); 282 cjmpL = &ctx->cjmpLR[0]; 283 cjmpR = &ctx->cjmpLR[dof]; 284 for (j=0; j<dof; j++) { 285 PetscScalar jmpL,jmpR; 286 jmpL = x[(i+0)*dof+j]-x[(i-1)*dof+j]; 287 jmpR = x[(i+1)*dof+j]-x[(i+0)*dof+j]; 288 for (k=0; k<dof; k++) { 289 cjmpL[k] += ctx->Rinv[k+j*dof]*jmpL; 290 cjmpR[k] += ctx->Rinv[k+j*dof]*jmpR; 291 } 292 } 293 /* Apply limiter to the left and right characteristic jumps */ 294 info.m = dof; 295 info.hx = hx; 296 (*ctx->limit)(&info,cjmpL,cjmpR,ctx->cslope); 297 for (j=0; j<dof; j++) ctx->cslope[j] /= hx; /* rescale to a slope */ 298 for (j=0; j<dof; j++) { 299 PetscScalar tmp = 0; 300 for (k=0; k<dof; k++) tmp += ctx->R[j+k*dof]*ctx->cslope[k]; 301 slope[i*dof+j] = tmp; 302 } 303 } 304 } 305 306 for (i=xs; i<xs+xm+1; i++) { 307 PetscReal maxspeed; 308 PetscScalar *uL,*uR; 309 if (i > len_slow) { 310 uL = &ctx->uLR[0]; 311 uR = &ctx->uLR[dof]; 312 for (j=0; j<dof; j++) { 313 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hx/2; 314 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hx/2; 315 } 316 CHKERRQ((*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed,ctx->xmin+hx*i,ctx->xmin,ctx->xmax)); 317 cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hx)); /* Max allowable value of 1/Delta t */ 318 if (i > xs) { 319 for (j=0; j<dof; j++) f[(i-len_slow-1)*dof+j] -= ctx->flux[j]/hx; 320 } 321 if (i < xs+xm) { 322 for (j=0; j<dof; j++) f[(i-len_slow)*dof+j] += ctx->flux[j]/hx; 323 } 324 } 325 if (i == len_slow) { 326 uL = &ctx->uLR[0]; 327 uR = &ctx->uLR[dof]; 328 for (j=0; j<dof; j++) { 329 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hx/2; 330 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hx/2; 331 } 332 CHKERRQ((*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed,ctx->xmin+hx*i,ctx->xmin,ctx->xmax)); 333 cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hx)); /* Max allowable value of 1/Delta t */ 334 if (i < xs+xm) { 335 for (j=0; j<dof; j++) f[(i-len_slow)*dof+j] += ctx->flux[j]/hx; 336 } 337 } 338 } 339 CHKERRQ(DMDAVecRestoreArray(da,Xloc,&x)); 340 CHKERRQ(VecRestoreArray(F,&f)); 341 CHKERRQ(DMDARestoreArray(da,PETSC_TRUE,&slope)); 342 CHKERRQ(DMRestoreLocalVector(da,&Xloc)); 343 PetscFunctionReturn(0); 344 } 345 346 PetscErrorCode FVRHSFunctionslow2(TS ts,PetscReal time,Vec X,Vec F,void *vctx) 347 { 348 FVCtx *ctx = (FVCtx*)vctx; 349 PetscInt i,j,k,Mx,dof,xs,xm,len_slow1,len_slow2; 350 PetscReal hx,cfl_idt = 0; 351 PetscScalar *x,*f,*slope; 352 Vec Xloc; 353 DM da; 354 355 PetscFunctionBeginUser; 356 CHKERRQ(TSGetDM(ts,&da)); 357 CHKERRQ(DMGetLocalVector(da,&Xloc)); 358 CHKERRQ(DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0)); 359 hx = (ctx->xmax-ctx->xmin)/Mx; 360 CHKERRQ(DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc)); 361 CHKERRQ(DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc)); 362 CHKERRQ(VecZeroEntries(F)); 363 CHKERRQ(DMDAVecGetArray(da,Xloc,&x)); 364 CHKERRQ(VecGetArray(F,&f)); 365 CHKERRQ(DMDAGetArray(da,PETSC_TRUE,&slope)); 366 CHKERRQ(DMDAGetCorners(da,&xs,0,0,&xm,0,0)); 367 CHKERRQ(ISGetSize(ctx->iss,&len_slow1)); 368 CHKERRQ(ISGetSize(ctx->iss2,&len_slow2)); 369 if (ctx->bctype == FVBC_OUTFLOW) { 370 for (i=xs-2; i<0; i++) { 371 for (j=0; j<dof; j++) x[i*dof+j] = x[j]; 372 } 373 for (i=Mx; i<xs+xm+2; i++) { 374 for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j]; 375 } 376 } 377 for (i=xs-1; i<xs+xm+1; i++) { 378 struct _LimitInfo info; 379 PetscScalar *cjmpL,*cjmpR; 380 if (i < len_slow1+len_slow2+1 && i > len_slow1-2) { 381 CHKERRQ((*ctx->physics.characteristic)(ctx->physics.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds,ctx->xmin+hx*i)); 382 CHKERRQ(PetscArrayzero(ctx->cjmpLR,2*dof)); 383 cjmpL = &ctx->cjmpLR[0]; 384 cjmpR = &ctx->cjmpLR[dof]; 385 for (j=0; j<dof; j++) { 386 PetscScalar jmpL,jmpR; 387 jmpL = x[(i+0)*dof+j]-x[(i-1)*dof+j]; 388 jmpR = x[(i+1)*dof+j]-x[(i+0)*dof+j]; 389 for (k=0; k<dof; k++) { 390 cjmpL[k] += ctx->Rinv[k+j*dof]*jmpL; 391 cjmpR[k] += ctx->Rinv[k+j*dof]*jmpR; 392 } 393 } 394 /* Apply limiter to the left and right characteristic jumps */ 395 info.m = dof; 396 info.hx = hx; 397 (*ctx->limit)(&info,cjmpL,cjmpR,ctx->cslope); 398 for (j=0; j<dof; j++) ctx->cslope[j] /= hx; /* rescale to a slope */ 399 for (j=0; j<dof; j++) { 400 PetscScalar tmp = 0; 401 for (k=0; k<dof; k++) tmp += ctx->R[j+k*dof]*ctx->cslope[k]; 402 slope[i*dof+j] = tmp; 403 } 404 } 405 } 406 407 for (i=xs; i<xs+xm+1; i++) { 408 PetscReal maxspeed; 409 PetscScalar *uL,*uR; 410 if (i < len_slow1+len_slow2 && i > len_slow1) { /* slow parts can be changed based on a */ 411 uL = &ctx->uLR[0]; 412 uR = &ctx->uLR[dof]; 413 for (j=0; j<dof; j++) { 414 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hx/2; 415 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hx/2; 416 } 417 CHKERRQ((*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed,ctx->xmin+hx*i,ctx->xmin,ctx->xmax)); 418 cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hx)); /* Max allowable value of 1/Delta t */ 419 if (i > xs) { 420 for (j=0; j<dof; j++) f[(i-len_slow1-1)*dof+j] -= ctx->flux[j]/hx; 421 } 422 if (i < xs+xm) { 423 for (j=0; j<dof; j++) f[(i-len_slow1)*dof+j] += ctx->flux[j]/hx; 424 } 425 } 426 if (i == len_slow1+len_slow2) { /* slow parts can be changed based on a */ 427 uL = &ctx->uLR[0]; 428 uR = &ctx->uLR[dof]; 429 for (j=0; j<dof; j++) { 430 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hx/2; 431 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hx/2; 432 } 433 CHKERRQ((*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed,ctx->xmin+hx*i,ctx->xmin,ctx->xmax)); 434 cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hx)); /* Max allowable value of 1/Delta t */ 435 if (i > xs) { 436 for (j=0; j<dof; j++) f[(i-len_slow1-1)*dof+j] -= ctx->flux[j]/hx; 437 } 438 } 439 if (i == len_slow1) { /* slow parts can be changed based on a */ 440 uL = &ctx->uLR[0]; 441 uR = &ctx->uLR[dof]; 442 for (j=0; j<dof; j++) { 443 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hx/2; 444 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hx/2; 445 } 446 CHKERRQ((*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed,ctx->xmin+hx*i,ctx->xmin,ctx->xmax)); 447 cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hx)); /* Max allowable value of 1/Delta t */ 448 if (i < xs+xm) { 449 for (j=0; j<dof; j++) f[(i-len_slow1)*dof+j] += ctx->flux[j]/hx; 450 } 451 } 452 } 453 454 CHKERRQ(DMDAVecRestoreArray(da,Xloc,&x)); 455 CHKERRQ(VecRestoreArray(F,&f)); 456 CHKERRQ(DMDARestoreArray(da,PETSC_TRUE,&slope)); 457 CHKERRQ(DMRestoreLocalVector(da,&Xloc)); 458 PetscFunctionReturn(0); 459 } 460 461 PetscErrorCode FVRHSFunctionfast2(TS ts,PetscReal time,Vec X,Vec F,void *vctx) 462 { 463 FVCtx *ctx = (FVCtx*)vctx; 464 PetscInt i,j,k,Mx,dof,xs,xm,len_slow1,len_slow2; 465 PetscReal hx,cfl_idt = 0; 466 PetscScalar *x,*f,*slope; 467 Vec Xloc; 468 DM da; 469 470 PetscFunctionBeginUser; 471 CHKERRQ(TSGetDM(ts,&da)); 472 CHKERRQ(DMGetLocalVector(da,&Xloc)); 473 CHKERRQ(DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0)); 474 hx = (ctx->xmax-ctx->xmin)/Mx; 475 CHKERRQ(DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc)); 476 CHKERRQ(DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc)); 477 CHKERRQ(VecZeroEntries(F)); 478 CHKERRQ(DMDAVecGetArray(da,Xloc,&x)); 479 CHKERRQ(VecGetArray(F,&f)); 480 CHKERRQ(DMDAGetArray(da,PETSC_TRUE,&slope)); 481 CHKERRQ(DMDAGetCorners(da,&xs,0,0,&xm,0,0)); 482 CHKERRQ(ISGetSize(ctx->iss,&len_slow1)); 483 CHKERRQ(ISGetSize(ctx->iss2,&len_slow2)); 484 485 if (ctx->bctype == FVBC_OUTFLOW) { 486 for (i=xs-2; i<0; i++) { 487 for (j=0; j<dof; j++) x[i*dof+j] = x[j]; 488 } 489 for (i=Mx; i<xs+xm+2; i++) { 490 for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j]; 491 } 492 } 493 for (i=xs-1; i<xs+xm+1; i++) { 494 struct _LimitInfo info; 495 PetscScalar *cjmpL,*cjmpR; 496 if (i > len_slow1+len_slow2-2) { 497 CHKERRQ((*ctx->physics.characteristic)(ctx->physics.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds,ctx->xmin+hx*i)); 498 CHKERRQ(PetscArrayzero(ctx->cjmpLR,2*dof)); 499 cjmpL = &ctx->cjmpLR[0]; 500 cjmpR = &ctx->cjmpLR[dof]; 501 for (j=0; j<dof; j++) { 502 PetscScalar jmpL,jmpR; 503 jmpL = x[(i+0)*dof+j]-x[(i-1)*dof+j]; 504 jmpR = x[(i+1)*dof+j]-x[(i+0)*dof+j]; 505 for (k=0; k<dof; k++) { 506 cjmpL[k] += ctx->Rinv[k+j*dof]*jmpL; 507 cjmpR[k] += ctx->Rinv[k+j*dof]*jmpR; 508 } 509 } 510 /* Apply limiter to the left and right characteristic jumps */ 511 info.m = dof; 512 info.hx = hx; 513 (*ctx->limit)(&info,cjmpL,cjmpR,ctx->cslope); 514 for (j=0; j<dof; j++) ctx->cslope[j] /= hx; /* rescale to a slope */ 515 for (j=0; j<dof; j++) { 516 PetscScalar tmp = 0; 517 for (k=0; k<dof; k++) tmp += ctx->R[j+k*dof]*ctx->cslope[k]; 518 slope[i*dof+j] = tmp; 519 } 520 } 521 } 522 523 for (i=xs; i<xs+xm+1; i++) { 524 PetscReal maxspeed; 525 PetscScalar *uL,*uR; 526 if (i > len_slow1+len_slow2) { 527 uL = &ctx->uLR[0]; 528 uR = &ctx->uLR[dof]; 529 for (j=0; j<dof; j++) { 530 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hx/2; 531 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hx/2; 532 } 533 CHKERRQ((*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed,ctx->xmin+hx*i,ctx->xmin,ctx->xmax)); 534 cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hx)); /* Max allowable value of 1/Delta t */ 535 if (i > xs) { 536 for (j=0; j<dof; j++) f[(i-len_slow1-len_slow2-1)*dof+j] -= ctx->flux[j]/hx; 537 } 538 if (i < xs+xm) { 539 for (j=0; j<dof; j++) f[(i-len_slow1-len_slow2)*dof+j] += ctx->flux[j]/hx; 540 } 541 } 542 if (i == len_slow1+len_slow2) { 543 uL = &ctx->uLR[0]; 544 uR = &ctx->uLR[dof]; 545 for (j=0; j<dof; j++) { 546 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hx/2; 547 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hx/2; 548 } 549 CHKERRQ((*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed,ctx->xmin+hx*i,ctx->xmin,ctx->xmax)); 550 cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hx)); /* Max allowable value of 1/Delta t */ 551 if (i < xs+xm) { 552 for (j=0; j<dof; j++) f[(i-len_slow1-len_slow2)*dof+j] += ctx->flux[j]/hx; 553 } 554 } 555 } 556 CHKERRQ(DMDAVecRestoreArray(da,Xloc,&x)); 557 CHKERRQ(VecRestoreArray(F,&f)); 558 CHKERRQ(DMDARestoreArray(da,PETSC_TRUE,&slope)); 559 CHKERRQ(DMRestoreLocalVector(da,&Xloc)); 560 PetscFunctionReturn(0); 561 } 562 563 int main(int argc,char *argv[]) 564 { 565 char lname[256] = "mc",physname[256] = "advect",final_fname[256] = "solution.m"; 566 PetscFunctionList limiters = 0,physics = 0; 567 MPI_Comm comm; 568 TS ts; 569 DM da; 570 Vec X,X0,R; 571 FVCtx ctx; 572 PetscInt i,k,dof,xs,xm,Mx,draw = 0,*index_slow,*index_fast,islow = 0,ifast = 0; 573 PetscBool view_final = PETSC_FALSE; 574 PetscReal ptime; 575 PetscErrorCode ierr; 576 577 CHKERRQ(PetscInitialize(&argc,&argv,0,help)); 578 comm = PETSC_COMM_WORLD; 579 CHKERRQ(PetscMemzero(&ctx,sizeof(ctx))); 580 581 /* Register limiters to be available on the command line */ 582 CHKERRQ(PetscFunctionListAdd(&limiters,"upwind" ,Limit_Upwind)); 583 CHKERRQ(PetscFunctionListAdd(&limiters,"lax-wendroff" ,Limit_LaxWendroff)); 584 CHKERRQ(PetscFunctionListAdd(&limiters,"beam-warming" ,Limit_BeamWarming)); 585 CHKERRQ(PetscFunctionListAdd(&limiters,"fromm" ,Limit_Fromm)); 586 CHKERRQ(PetscFunctionListAdd(&limiters,"minmod" ,Limit_Minmod)); 587 CHKERRQ(PetscFunctionListAdd(&limiters,"superbee" ,Limit_Superbee)); 588 CHKERRQ(PetscFunctionListAdd(&limiters,"mc" ,Limit_MC)); 589 CHKERRQ(PetscFunctionListAdd(&limiters,"vanleer" ,Limit_VanLeer)); 590 CHKERRQ(PetscFunctionListAdd(&limiters,"vanalbada" ,Limit_VanAlbada)); 591 CHKERRQ(PetscFunctionListAdd(&limiters,"vanalbadatvd" ,Limit_VanAlbadaTVD)); 592 CHKERRQ(PetscFunctionListAdd(&limiters,"koren" ,Limit_Koren)); 593 CHKERRQ(PetscFunctionListAdd(&limiters,"korensym" ,Limit_KorenSym)); 594 CHKERRQ(PetscFunctionListAdd(&limiters,"koren3" ,Limit_Koren3)); 595 CHKERRQ(PetscFunctionListAdd(&limiters,"cada-torrilhon2" ,Limit_CadaTorrilhon2)); 596 CHKERRQ(PetscFunctionListAdd(&limiters,"cada-torrilhon3-r0p1",Limit_CadaTorrilhon3R0p1)); 597 CHKERRQ(PetscFunctionListAdd(&limiters,"cada-torrilhon3-r1" ,Limit_CadaTorrilhon3R1)); 598 CHKERRQ(PetscFunctionListAdd(&limiters,"cada-torrilhon3-r10" ,Limit_CadaTorrilhon3R10)); 599 CHKERRQ(PetscFunctionListAdd(&limiters,"cada-torrilhon3-r100",Limit_CadaTorrilhon3R100)); 600 601 /* Register physical models to be available on the command line */ 602 CHKERRQ(PetscFunctionListAdd(&physics,"advect" ,PhysicsCreate_Advect)); 603 604 ctx.comm = comm; 605 ctx.cfl = 0.9; 606 ctx.bctype = FVBC_PERIODIC; 607 ctx.xmin = -1.0; 608 ctx.xmax = 1.0; 609 ierr = PetscOptionsBegin(comm,NULL,"Finite Volume solver options","");CHKERRQ(ierr); 610 CHKERRQ(PetscOptionsReal("-xmin","X min","",ctx.xmin,&ctx.xmin,NULL)); 611 CHKERRQ(PetscOptionsReal("-xmax","X max","",ctx.xmax,&ctx.xmax,NULL)); 612 CHKERRQ(PetscOptionsFList("-limit","Name of flux limiter to use","",limiters,lname,lname,sizeof(lname),NULL)); 613 CHKERRQ(PetscOptionsInt("-draw","Draw solution vector, bitwise OR of (1=initial,2=final,4=final error)","",draw,&draw,NULL)); 614 CHKERRQ(PetscOptionsString("-view_final","Write final solution in ASCII MATLAB format to given file name","",final_fname,final_fname,sizeof(final_fname),&view_final)); 615 CHKERRQ(PetscOptionsInt("-initial","Initial condition (depends on the physics)","",ctx.initial,&ctx.initial,NULL)); 616 CHKERRQ(PetscOptionsBool("-simulation","Compare errors with reference solution","",ctx.simulation,&ctx.simulation,NULL)); 617 CHKERRQ(PetscOptionsReal("-cfl","CFL number to time step at","",ctx.cfl,&ctx.cfl,NULL)); 618 CHKERRQ(PetscOptionsEnum("-bc_type","Boundary condition","",FVBCTypes,(PetscEnum)ctx.bctype,(PetscEnum*)&ctx.bctype,NULL)); 619 CHKERRQ(PetscOptionsBool("-recursive_split","Split the domain recursively","",ctx.recursive,&ctx.recursive,NULL)); 620 ierr = PetscOptionsEnd();CHKERRQ(ierr); 621 622 /* Choose the limiter from the list of registered limiters */ 623 CHKERRQ(PetscFunctionListFind(limiters,lname,&ctx.limit)); 624 PetscCheck(ctx.limit,PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Limiter '%s' not found",lname); 625 626 /* Choose the physics from the list of registered models */ 627 { 628 PetscErrorCode (*r)(FVCtx*); 629 CHKERRQ(PetscFunctionListFind(physics,physname,&r)); 630 PetscCheck(r,PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Physics '%s' not found",physname); 631 /* Create the physics, will set the number of fields and their names */ 632 CHKERRQ((*r)(&ctx)); 633 } 634 635 /* Create a DMDA to manage the parallel grid */ 636 CHKERRQ(DMDACreate1d(comm,DM_BOUNDARY_PERIODIC,50,ctx.physics.dof,2,NULL,&da)); 637 CHKERRQ(DMSetFromOptions(da)); 638 CHKERRQ(DMSetUp(da)); 639 /* Inform the DMDA of the field names provided by the physics. */ 640 /* The names will be shown in the title bars when run with -ts_monitor_draw_solution */ 641 for (i=0; i<ctx.physics.dof; i++) { 642 CHKERRQ(DMDASetFieldName(da,i,ctx.physics.fieldname[i])); 643 } 644 CHKERRQ(DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0)); 645 CHKERRQ(DMDAGetCorners(da,&xs,0,0,&xm,0,0)); 646 647 /* Set coordinates of cell centers */ 648 CHKERRQ(DMDASetUniformCoordinates(da,ctx.xmin+0.5*(ctx.xmax-ctx.xmin)/Mx,ctx.xmax+0.5*(ctx.xmax-ctx.xmin)/Mx,0,0,0,0)); 649 650 /* Allocate work space for the Finite Volume solver (so it doesn't have to be reallocated on each function evaluation) */ 651 CHKERRQ(PetscMalloc4(dof*dof,&ctx.R,dof*dof,&ctx.Rinv,2*dof,&ctx.cjmpLR,1*dof,&ctx.cslope)); 652 CHKERRQ(PetscMalloc3(2*dof,&ctx.uLR,dof,&ctx.flux,dof,&ctx.speeds)); 653 654 /* Create a vector to store the solution and to save the initial state */ 655 CHKERRQ(DMCreateGlobalVector(da,&X)); 656 CHKERRQ(VecDuplicate(X,&X0)); 657 CHKERRQ(VecDuplicate(X,&R)); 658 659 /*-------------------------------- create index for slow parts and fast parts ----------------------------------------*/ 660 CHKERRQ(PetscMalloc1(xm*dof,&index_slow)); 661 CHKERRQ(PetscMalloc1(xm*dof,&index_fast)); 662 for (i=xs; i<xs+xm; i++) { 663 if (ctx.xmin+i*(ctx.xmax-ctx.xmin)/(PetscReal)Mx+0.5*(ctx.xmax-ctx.xmin)/(PetscReal)Mx < 0) 664 for (k=0; k<dof; k++) index_slow[islow++] = i*dof+k; 665 else 666 for (k=0; k<dof; k++) index_fast[ifast++] = i*dof+k; 667 } /* this step need to be changed based on discontinuous point of a */ 668 CHKERRQ(ISCreateGeneral(PETSC_COMM_WORLD,islow,index_slow,PETSC_COPY_VALUES,&ctx.iss)); 669 CHKERRQ(ISCreateGeneral(PETSC_COMM_WORLD,ifast,index_fast,PETSC_COPY_VALUES,&ctx.isf)); 670 671 /* Create a time-stepping object */ 672 CHKERRQ(TSCreate(comm,&ts)); 673 CHKERRQ(TSSetDM(ts,da)); 674 CHKERRQ(TSSetRHSFunction(ts,R,FVRHSFunction,&ctx)); 675 CHKERRQ(TSRHSSplitSetIS(ts,"slow",ctx.iss)); 676 CHKERRQ(TSRHSSplitSetIS(ts,"fast",ctx.isf)); 677 CHKERRQ(TSRHSSplitSetRHSFunction(ts,"slow",NULL,FVRHSFunctionslow,&ctx)); 678 CHKERRQ(TSRHSSplitSetRHSFunction(ts,"fast",NULL,FVRHSFunctionfast,&ctx)); 679 680 if (ctx.recursive) { 681 TS subts; 682 islow = 0; 683 ifast = 0; 684 for (i=xs; i<xs+xm; i++) { 685 PetscReal coord = ctx.xmin+i*(ctx.xmax-ctx.xmin)/(PetscReal)Mx+0.5*(ctx.xmax-ctx.xmin)/(PetscReal)Mx; 686 if (coord >= 0 && coord < ctx.xmin+(ctx.xmax-ctx.xmin)*3/4.) 687 for (k=0; k<dof; k++) index_slow[islow++] = i*dof+k; 688 if (coord >= ctx.xmin+(ctx.xmax-ctx.xmin)*3/4.) 689 for (k=0; k<dof; k++) index_fast[ifast++] = i*dof+k; 690 } /* this step need to be changed based on discontinuous point of a */ 691 CHKERRQ(ISCreateGeneral(PETSC_COMM_WORLD,islow,index_slow,PETSC_COPY_VALUES,&ctx.iss2)); 692 CHKERRQ(ISCreateGeneral(PETSC_COMM_WORLD,ifast,index_fast,PETSC_COPY_VALUES,&ctx.isf2)); 693 694 CHKERRQ(TSRHSSplitGetSubTS(ts,"fast",&subts)); 695 CHKERRQ(TSRHSSplitSetIS(subts,"slow",ctx.iss2)); 696 CHKERRQ(TSRHSSplitSetIS(subts,"fast",ctx.isf2)); 697 CHKERRQ(TSRHSSplitSetRHSFunction(subts,"slow",NULL,FVRHSFunctionslow2,&ctx)); 698 CHKERRQ(TSRHSSplitSetRHSFunction(subts,"fast",NULL,FVRHSFunctionfast2,&ctx)); 699 } 700 701 /*CHKERRQ(TSSetType(ts,TSSSP));*/ 702 CHKERRQ(TSSetType(ts,TSMPRK)); 703 CHKERRQ(TSSetMaxTime(ts,10)); 704 CHKERRQ(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER)); 705 706 /* Compute initial conditions and starting time step */ 707 CHKERRQ(FVSample(&ctx,da,0,X0)); 708 CHKERRQ(FVRHSFunction(ts,0,X0,X,(void*)&ctx)); /* Initial function evaluation, only used to determine max speed */ 709 CHKERRQ(VecCopy(X0,X)); /* The function value was not used so we set X=X0 again */ 710 CHKERRQ(TSSetTimeStep(ts,ctx.cfl/ctx.cfl_idt)); 711 CHKERRQ(TSSetFromOptions(ts)); /* Take runtime options */ 712 CHKERRQ(SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD)); 713 { 714 PetscInt steps; 715 PetscScalar mass_initial,mass_final,mass_difference; 716 717 CHKERRQ(TSSolve(ts,X)); 718 CHKERRQ(TSGetSolveTime(ts,&ptime)); 719 CHKERRQ(TSGetStepNumber(ts,&steps)); 720 CHKERRQ(PetscPrintf(comm,"Final time %g, steps %D\n",(double)ptime,steps)); 721 /* calculate the total mass at initial time and final time */ 722 mass_initial = 0.0; 723 mass_final = 0.0; 724 CHKERRQ(VecSum(X0,&mass_initial)); 725 CHKERRQ(VecSum(X,&mass_final)); 726 mass_difference = (ctx.xmax-ctx.xmin)/(PetscScalar)Mx*(mass_final - mass_initial); 727 CHKERRQ(PetscPrintf(comm,"Mass difference %g\n",(double)mass_difference)); 728 if (ctx.simulation) { 729 PetscViewer fd; 730 char filename[PETSC_MAX_PATH_LEN] = "binaryoutput"; 731 Vec XR; 732 PetscReal nrm1,nrmsup; 733 PetscBool flg; 734 735 CHKERRQ(PetscOptionsGetString(NULL,NULL,"-f",filename,sizeof(filename),&flg)); 736 PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_USER,"Must indicate binary file with the -f option"); 737 CHKERRQ(PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&fd)); 738 CHKERRQ(VecDuplicate(X0,&XR)); 739 CHKERRQ(VecLoad(XR,fd)); 740 CHKERRQ(PetscViewerDestroy(&fd)); 741 CHKERRQ(VecAYPX(XR,-1,X)); 742 CHKERRQ(VecNorm(XR,NORM_1,&nrm1)); 743 CHKERRQ(VecNorm(XR,NORM_INFINITY,&nrmsup)); 744 nrm1 /= Mx; 745 CHKERRQ(PetscPrintf(comm,"Error ||x-x_e||_1 %g ||x-x_e||_sup %g\n",(double)nrm1,(double)nrmsup)); 746 CHKERRQ(VecDestroy(&XR)); 747 } 748 } 749 750 CHKERRQ(SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD)); 751 if (draw & 0x1) CHKERRQ(VecView(X0,PETSC_VIEWER_DRAW_WORLD)); 752 if (draw & 0x2) CHKERRQ(VecView(X,PETSC_VIEWER_DRAW_WORLD)); 753 if (draw & 0x4) { 754 Vec Y; 755 CHKERRQ(VecDuplicate(X,&Y)); 756 CHKERRQ(FVSample(&ctx,da,ptime,Y)); 757 CHKERRQ(VecAYPX(Y,-1,X)); 758 CHKERRQ(VecView(Y,PETSC_VIEWER_DRAW_WORLD)); 759 CHKERRQ(VecDestroy(&Y)); 760 } 761 762 if (view_final) { 763 PetscViewer viewer; 764 CHKERRQ(PetscViewerASCIIOpen(PETSC_COMM_WORLD,final_fname,&viewer)); 765 CHKERRQ(PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB)); 766 CHKERRQ(VecView(X,viewer)); 767 CHKERRQ(PetscViewerPopFormat(viewer)); 768 CHKERRQ(PetscViewerDestroy(&viewer)); 769 } 770 771 /* Clean up */ 772 CHKERRQ((*ctx.physics.destroy)(ctx.physics.user)); 773 for (i=0; i<ctx.physics.dof; i++) CHKERRQ(PetscFree(ctx.physics.fieldname[i])); 774 CHKERRQ(PetscFree4(ctx.R,ctx.Rinv,ctx.cjmpLR,ctx.cslope)); 775 CHKERRQ(PetscFree3(ctx.uLR,ctx.flux,ctx.speeds)); 776 CHKERRQ(ISDestroy(&ctx.iss)); 777 CHKERRQ(ISDestroy(&ctx.isf)); 778 CHKERRQ(ISDestroy(&ctx.iss2)); 779 CHKERRQ(ISDestroy(&ctx.isf2)); 780 CHKERRQ(VecDestroy(&X)); 781 CHKERRQ(VecDestroy(&X0)); 782 CHKERRQ(VecDestroy(&R)); 783 CHKERRQ(DMDestroy(&da)); 784 CHKERRQ(TSDestroy(&ts)); 785 CHKERRQ(PetscFree(index_slow)); 786 CHKERRQ(PetscFree(index_fast)); 787 CHKERRQ(PetscFunctionListDestroy(&limiters)); 788 CHKERRQ(PetscFunctionListDestroy(&physics)); 789 CHKERRQ(PetscFinalize()); 790 return 0; 791 } 792 793 /*TEST 794 build: 795 requires: !complex 796 depends: finitevolume1d.c 797 798 test: 799 args: -da_grid_x 60 -initial 1 -xmin -1 -xmax 1 -limit mc -physics_advect_a1 1 -physics_advect_a2 2 -ts_dt 0.025 -ts_max_steps 24 -ts_type rk -ts_rk_type 2a -ts_rk_dtratio 2 -ts_rk_multirate -ts_use_splitrhsfunction 0 800 801 test: 802 suffix: 2 803 args: -da_grid_x 60 -initial 1 -xmin -1 -xmax 1 -limit mc -physics_advect_a1 1 -physics_advect_a2 2 -ts_dt 0.025 -ts_max_steps 24 -ts_type rk -ts_rk_type 2a -ts_rk_dtratio 2 -ts_rk_multirate -ts_use_splitrhsfunction 1 804 output_file: output/ex5_1.out 805 806 test: 807 suffix: 3 808 args: -da_grid_x 60 -initial 1 -xmin -1 -xmax 1 -limit mc -physics_advect_a1 1 -physics_advect_a2 2 -ts_dt 0.025 -ts_max_steps 24 -ts_type mprk -ts_mprk_type 2a22 -ts_use_splitrhsfunction 0 809 810 test: 811 suffix: 4 812 args: -da_grid_x 60 -initial 1 -xmin -1 -xmax 1 -limit mc -physics_advect_a1 1 -physics_advect_a2 2 -ts_dt 0.025 -ts_max_steps 24 -ts_type mprk -ts_mprk_type 2a22 -ts_use_splitrhsfunction 1 813 output_file: output/ex5_3.out 814 815 test: 816 suffix: 5 817 args: -da_grid_x 60 -initial 1 -xmin -1 -xmax 1 -limit mc -physics_advect_a1 1 -physics_advect_a2 2 -ts_dt 0.025 -ts_max_steps 24 -ts_type rk -ts_rk_type 2a -ts_rk_dtratio 2 -ts_rk_multirate -ts_use_splitrhsfunction 0 -recursive_split 818 819 test: 820 suffix: 6 821 args: -da_grid_x 60 -initial 1 -xmin -1 -xmax 1 -limit mc -physics_advect_a1 1 -physics_advect_a2 2 -ts_dt 0.025 -ts_max_steps 24 -ts_type rk -ts_rk_type 2a -ts_rk_dtratio 2 -ts_rk_multirate -ts_use_splitrhsfunction 1 -recursive_split 822 TEST*/ 823