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