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 PetscFunctionBeginUser; 576 PetscCall(PetscInitialize(&argc,&argv,0,help)); 577 comm = PETSC_COMM_WORLD; 578 PetscCall(PetscMemzero(&ctx,sizeof(ctx))); 579 580 /* Register limiters to be available on the command line */ 581 PetscCall(PetscFunctionListAdd(&limiters,"upwind" ,Limit_Upwind)); 582 PetscCall(PetscFunctionListAdd(&limiters,"lax-wendroff" ,Limit_LaxWendroff)); 583 PetscCall(PetscFunctionListAdd(&limiters,"beam-warming" ,Limit_BeamWarming)); 584 PetscCall(PetscFunctionListAdd(&limiters,"fromm" ,Limit_Fromm)); 585 PetscCall(PetscFunctionListAdd(&limiters,"minmod" ,Limit_Minmod)); 586 PetscCall(PetscFunctionListAdd(&limiters,"superbee" ,Limit_Superbee)); 587 PetscCall(PetscFunctionListAdd(&limiters,"mc" ,Limit_MC)); 588 PetscCall(PetscFunctionListAdd(&limiters,"vanleer" ,Limit_VanLeer)); 589 PetscCall(PetscFunctionListAdd(&limiters,"vanalbada" ,Limit_VanAlbada)); 590 PetscCall(PetscFunctionListAdd(&limiters,"vanalbadatvd" ,Limit_VanAlbadaTVD)); 591 PetscCall(PetscFunctionListAdd(&limiters,"koren" ,Limit_Koren)); 592 PetscCall(PetscFunctionListAdd(&limiters,"korensym" ,Limit_KorenSym)); 593 PetscCall(PetscFunctionListAdd(&limiters,"koren3" ,Limit_Koren3)); 594 PetscCall(PetscFunctionListAdd(&limiters,"cada-torrilhon2" ,Limit_CadaTorrilhon2)); 595 PetscCall(PetscFunctionListAdd(&limiters,"cada-torrilhon3-r0p1",Limit_CadaTorrilhon3R0p1)); 596 PetscCall(PetscFunctionListAdd(&limiters,"cada-torrilhon3-r1" ,Limit_CadaTorrilhon3R1)); 597 PetscCall(PetscFunctionListAdd(&limiters,"cada-torrilhon3-r10" ,Limit_CadaTorrilhon3R10)); 598 PetscCall(PetscFunctionListAdd(&limiters,"cada-torrilhon3-r100",Limit_CadaTorrilhon3R100)); 599 600 /* Register physical models to be available on the command line */ 601 PetscCall(PetscFunctionListAdd(&physics,"advect" ,PhysicsCreate_Advect)); 602 603 ctx.comm = comm; 604 ctx.cfl = 0.9; 605 ctx.bctype = FVBC_PERIODIC; 606 ctx.xmin = -1.0; 607 ctx.xmax = 1.0; 608 PetscOptionsBegin(comm,NULL,"Finite Volume solver options",""); 609 PetscCall(PetscOptionsReal("-xmin","X min","",ctx.xmin,&ctx.xmin,NULL)); 610 PetscCall(PetscOptionsReal("-xmax","X max","",ctx.xmax,&ctx.xmax,NULL)); 611 PetscCall(PetscOptionsFList("-limit","Name of flux limiter to use","",limiters,lname,lname,sizeof(lname),NULL)); 612 PetscCall(PetscOptionsInt("-draw","Draw solution vector, bitwise OR of (1=initial,2=final,4=final error)","",draw,&draw,NULL)); 613 PetscCall(PetscOptionsString("-view_final","Write final solution in ASCII MATLAB format to given file name","",final_fname,final_fname,sizeof(final_fname),&view_final)); 614 PetscCall(PetscOptionsInt("-initial","Initial condition (depends on the physics)","",ctx.initial,&ctx.initial,NULL)); 615 PetscCall(PetscOptionsBool("-simulation","Compare errors with reference solution","",ctx.simulation,&ctx.simulation,NULL)); 616 PetscCall(PetscOptionsReal("-cfl","CFL number to time step at","",ctx.cfl,&ctx.cfl,NULL)); 617 PetscCall(PetscOptionsEnum("-bc_type","Boundary condition","",FVBCTypes,(PetscEnum)ctx.bctype,(PetscEnum*)&ctx.bctype,NULL)); 618 PetscCall(PetscOptionsBool("-recursive_split","Split the domain recursively","",ctx.recursive,&ctx.recursive,NULL)); 619 PetscOptionsEnd(); 620 621 /* Choose the limiter from the list of registered limiters */ 622 PetscCall(PetscFunctionListFind(limiters,lname,&ctx.limit)); 623 PetscCheck(ctx.limit,PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Limiter '%s' not found",lname); 624 625 /* Choose the physics from the list of registered models */ 626 { 627 PetscErrorCode (*r)(FVCtx*); 628 PetscCall(PetscFunctionListFind(physics,physname,&r)); 629 PetscCheck(r,PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Physics '%s' not found",physname); 630 /* Create the physics, will set the number of fields and their names */ 631 PetscCall((*r)(&ctx)); 632 } 633 634 /* Create a DMDA to manage the parallel grid */ 635 PetscCall(DMDACreate1d(comm,DM_BOUNDARY_PERIODIC,50,ctx.physics.dof,2,NULL,&da)); 636 PetscCall(DMSetFromOptions(da)); 637 PetscCall(DMSetUp(da)); 638 /* Inform the DMDA of the field names provided by the physics. */ 639 /* The names will be shown in the title bars when run with -ts_monitor_draw_solution */ 640 for (i=0; i<ctx.physics.dof; i++) { 641 PetscCall(DMDASetFieldName(da,i,ctx.physics.fieldname[i])); 642 } 643 PetscCall(DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0)); 644 PetscCall(DMDAGetCorners(da,&xs,0,0,&xm,0,0)); 645 646 /* Set coordinates of cell centers */ 647 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)); 648 649 /* Allocate work space for the Finite Volume solver (so it doesn't have to be reallocated on each function evaluation) */ 650 PetscCall(PetscMalloc4(dof*dof,&ctx.R,dof*dof,&ctx.Rinv,2*dof,&ctx.cjmpLR,1*dof,&ctx.cslope)); 651 PetscCall(PetscMalloc3(2*dof,&ctx.uLR,dof,&ctx.flux,dof,&ctx.speeds)); 652 653 /* Create a vector to store the solution and to save the initial state */ 654 PetscCall(DMCreateGlobalVector(da,&X)); 655 PetscCall(VecDuplicate(X,&X0)); 656 PetscCall(VecDuplicate(X,&R)); 657 658 /*-------------------------------- create index for slow parts and fast parts ----------------------------------------*/ 659 PetscCall(PetscMalloc1(xm*dof,&index_slow)); 660 PetscCall(PetscMalloc1(xm*dof,&index_fast)); 661 for (i=xs; i<xs+xm; i++) { 662 if (ctx.xmin+i*(ctx.xmax-ctx.xmin)/(PetscReal)Mx+0.5*(ctx.xmax-ctx.xmin)/(PetscReal)Mx < 0) 663 for (k=0; k<dof; k++) index_slow[islow++] = i*dof+k; 664 else 665 for (k=0; k<dof; k++) index_fast[ifast++] = i*dof+k; 666 } /* this step need to be changed based on discontinuous point of a */ 667 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD,islow,index_slow,PETSC_COPY_VALUES,&ctx.iss)); 668 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD,ifast,index_fast,PETSC_COPY_VALUES,&ctx.isf)); 669 670 /* Create a time-stepping object */ 671 PetscCall(TSCreate(comm,&ts)); 672 PetscCall(TSSetDM(ts,da)); 673 PetscCall(TSSetRHSFunction(ts,R,FVRHSFunction,&ctx)); 674 PetscCall(TSRHSSplitSetIS(ts,"slow",ctx.iss)); 675 PetscCall(TSRHSSplitSetIS(ts,"fast",ctx.isf)); 676 PetscCall(TSRHSSplitSetRHSFunction(ts,"slow",NULL,FVRHSFunctionslow,&ctx)); 677 PetscCall(TSRHSSplitSetRHSFunction(ts,"fast",NULL,FVRHSFunctionfast,&ctx)); 678 679 if (ctx.recursive) { 680 TS subts; 681 islow = 0; 682 ifast = 0; 683 for (i=xs; i<xs+xm; i++) { 684 PetscReal coord = ctx.xmin+i*(ctx.xmax-ctx.xmin)/(PetscReal)Mx+0.5*(ctx.xmax-ctx.xmin)/(PetscReal)Mx; 685 if (coord >= 0 && coord < ctx.xmin+(ctx.xmax-ctx.xmin)*3/4.) 686 for (k=0; k<dof; k++) index_slow[islow++] = i*dof+k; 687 if (coord >= ctx.xmin+(ctx.xmax-ctx.xmin)*3/4.) 688 for (k=0; k<dof; k++) index_fast[ifast++] = i*dof+k; 689 } /* this step need to be changed based on discontinuous point of a */ 690 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD,islow,index_slow,PETSC_COPY_VALUES,&ctx.iss2)); 691 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD,ifast,index_fast,PETSC_COPY_VALUES,&ctx.isf2)); 692 693 PetscCall(TSRHSSplitGetSubTS(ts,"fast",&subts)); 694 PetscCall(TSRHSSplitSetIS(subts,"slow",ctx.iss2)); 695 PetscCall(TSRHSSplitSetIS(subts,"fast",ctx.isf2)); 696 PetscCall(TSRHSSplitSetRHSFunction(subts,"slow",NULL,FVRHSFunctionslow2,&ctx)); 697 PetscCall(TSRHSSplitSetRHSFunction(subts,"fast",NULL,FVRHSFunctionfast2,&ctx)); 698 } 699 700 /*PetscCall(TSSetType(ts,TSSSP));*/ 701 PetscCall(TSSetType(ts,TSMPRK)); 702 PetscCall(TSSetMaxTime(ts,10)); 703 PetscCall(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER)); 704 705 /* Compute initial conditions and starting time step */ 706 PetscCall(FVSample(&ctx,da,0,X0)); 707 PetscCall(FVRHSFunction(ts,0,X0,X,(void*)&ctx)); /* Initial function evaluation, only used to determine max speed */ 708 PetscCall(VecCopy(X0,X)); /* The function value was not used so we set X=X0 again */ 709 PetscCall(TSSetTimeStep(ts,ctx.cfl/ctx.cfl_idt)); 710 PetscCall(TSSetFromOptions(ts)); /* Take runtime options */ 711 PetscCall(SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD)); 712 { 713 PetscInt steps; 714 PetscScalar mass_initial,mass_final,mass_difference; 715 716 PetscCall(TSSolve(ts,X)); 717 PetscCall(TSGetSolveTime(ts,&ptime)); 718 PetscCall(TSGetStepNumber(ts,&steps)); 719 PetscCall(PetscPrintf(comm,"Final time %g, steps %" PetscInt_FMT "\n",(double)ptime,steps)); 720 /* calculate the total mass at initial time and final time */ 721 mass_initial = 0.0; 722 mass_final = 0.0; 723 PetscCall(VecSum(X0,&mass_initial)); 724 PetscCall(VecSum(X,&mass_final)); 725 mass_difference = (ctx.xmax-ctx.xmin)/(PetscScalar)Mx*(mass_final - mass_initial); 726 PetscCall(PetscPrintf(comm,"Mass difference %g\n",(double)mass_difference)); 727 if (ctx.simulation) { 728 PetscViewer fd; 729 char filename[PETSC_MAX_PATH_LEN] = "binaryoutput"; 730 Vec XR; 731 PetscReal nrm1,nrmsup; 732 PetscBool flg; 733 734 PetscCall(PetscOptionsGetString(NULL,NULL,"-f",filename,sizeof(filename),&flg)); 735 PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_USER,"Must indicate binary file with the -f option"); 736 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&fd)); 737 PetscCall(VecDuplicate(X0,&XR)); 738 PetscCall(VecLoad(XR,fd)); 739 PetscCall(PetscViewerDestroy(&fd)); 740 PetscCall(VecAYPX(XR,-1,X)); 741 PetscCall(VecNorm(XR,NORM_1,&nrm1)); 742 PetscCall(VecNorm(XR,NORM_INFINITY,&nrmsup)); 743 nrm1 /= Mx; 744 PetscCall(PetscPrintf(comm,"Error ||x-x_e||_1 %g ||x-x_e||_sup %g\n",(double)nrm1,(double)nrmsup)); 745 PetscCall(VecDestroy(&XR)); 746 } 747 } 748 749 PetscCall(SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD)); 750 if (draw & 0x1) PetscCall(VecView(X0,PETSC_VIEWER_DRAW_WORLD)); 751 if (draw & 0x2) PetscCall(VecView(X,PETSC_VIEWER_DRAW_WORLD)); 752 if (draw & 0x4) { 753 Vec Y; 754 PetscCall(VecDuplicate(X,&Y)); 755 PetscCall(FVSample(&ctx,da,ptime,Y)); 756 PetscCall(VecAYPX(Y,-1,X)); 757 PetscCall(VecView(Y,PETSC_VIEWER_DRAW_WORLD)); 758 PetscCall(VecDestroy(&Y)); 759 } 760 761 if (view_final) { 762 PetscViewer viewer; 763 PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD,final_fname,&viewer)); 764 PetscCall(PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB)); 765 PetscCall(VecView(X,viewer)); 766 PetscCall(PetscViewerPopFormat(viewer)); 767 PetscCall(PetscViewerDestroy(&viewer)); 768 } 769 770 /* Clean up */ 771 PetscCall((*ctx.physics.destroy)(ctx.physics.user)); 772 for (i=0; i<ctx.physics.dof; i++) PetscCall(PetscFree(ctx.physics.fieldname[i])); 773 PetscCall(PetscFree4(ctx.R,ctx.Rinv,ctx.cjmpLR,ctx.cslope)); 774 PetscCall(PetscFree3(ctx.uLR,ctx.flux,ctx.speeds)); 775 PetscCall(ISDestroy(&ctx.iss)); 776 PetscCall(ISDestroy(&ctx.isf)); 777 PetscCall(ISDestroy(&ctx.iss2)); 778 PetscCall(ISDestroy(&ctx.isf2)); 779 PetscCall(VecDestroy(&X)); 780 PetscCall(VecDestroy(&X0)); 781 PetscCall(VecDestroy(&R)); 782 PetscCall(DMDestroy(&da)); 783 PetscCall(TSDestroy(&ts)); 784 PetscCall(PetscFree(index_slow)); 785 PetscCall(PetscFree(index_fast)); 786 PetscCall(PetscFunctionListDestroy(&limiters)); 787 PetscCall(PetscFunctionListDestroy(&physics)); 788 PetscCall(PetscFinalize()); 789 return 0; 790 } 791 792 /*TEST 793 build: 794 requires: !complex 795 depends: finitevolume1d.c 796 797 test: 798 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 799 800 test: 801 suffix: 2 802 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 803 output_file: output/ex5_1.out 804 805 test: 806 suffix: 3 807 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 808 809 test: 810 suffix: 4 811 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 812 output_file: output/ex5_3.out 813 814 test: 815 suffix: 5 816 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 817 818 test: 819 suffix: 6 820 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 821 TEST*/ 822