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 PetscCall(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 PetscCall(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","");PetscCall(ierr); 130 { 131 PetscCall(PetscOptionsReal("-physics_advect_a1","Speed1","",user->a[0],&user->a[0],NULL)); 132 PetscCall(PetscOptionsReal("-physics_advect_a2","Speed2","",user->a[1],&user->a[1],NULL)); 133 } 134 ierr = PetscOptionsEnd();PetscCall(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 PetscCall(TSGetDM(ts,&da)); 151 PetscCall(DMGetLocalVector(da,&Xloc)); 152 PetscCall(DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0)); 153 hx = (ctx->xmax-ctx->xmin)/Mx; 154 PetscCall(DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc)); 155 PetscCall(DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc)); 156 PetscCall(VecZeroEntries(F)); 157 PetscCall(DMDAVecGetArray(da,Xloc,&x)); 158 PetscCall(VecGetArray(F,&f)); 159 PetscCall(DMDAGetArray(da,PETSC_TRUE,&slope)); 160 PetscCall(DMDAGetCorners(da,&xs,0,0,&xm,0,0)); 161 PetscCall(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 PetscCall((*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 PetscCall(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 PetscCall((*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 PetscCall((*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 PetscCall(DMDAVecRestoreArray(da,Xloc,&x)); 237 PetscCall(VecRestoreArray(F,&f)); 238 PetscCall(DMDARestoreArray(da,PETSC_TRUE,&slope)); 239 PetscCall(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 PetscCall(TSGetDM(ts,&da)); 256 PetscCall(DMGetLocalVector(da,&Xloc)); 257 PetscCall(DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0)); 258 hx = (ctx->xmax-ctx->xmin)/Mx; 259 PetscCall(DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc)); 260 PetscCall(DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc)); 261 PetscCall(VecZeroEntries(F)); 262 PetscCall(DMDAVecGetArray(da,Xloc,&x)); 263 PetscCall(VecGetArray(F,&f)); 264 PetscCall(DMDAGetArray(da,PETSC_TRUE,&slope)); 265 PetscCall(DMDAGetCorners(da,&xs,0,0,&xm,0,0)); 266 PetscCall(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 PetscCall((*ctx->physics.characteristic)(ctx->physics.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds,ctx->xmin+hx*i)); 281 PetscCall(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 PetscCall((*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 PetscCall((*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 PetscCall(DMDAVecRestoreArray(da,Xloc,&x)); 340 PetscCall(VecRestoreArray(F,&f)); 341 PetscCall(DMDARestoreArray(da,PETSC_TRUE,&slope)); 342 PetscCall(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 PetscCall(TSGetDM(ts,&da)); 357 PetscCall(DMGetLocalVector(da,&Xloc)); 358 PetscCall(DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0)); 359 hx = (ctx->xmax-ctx->xmin)/Mx; 360 PetscCall(DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc)); 361 PetscCall(DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc)); 362 PetscCall(VecZeroEntries(F)); 363 PetscCall(DMDAVecGetArray(da,Xloc,&x)); 364 PetscCall(VecGetArray(F,&f)); 365 PetscCall(DMDAGetArray(da,PETSC_TRUE,&slope)); 366 PetscCall(DMDAGetCorners(da,&xs,0,0,&xm,0,0)); 367 PetscCall(ISGetSize(ctx->iss,&len_slow1)); 368 PetscCall(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 PetscCall((*ctx->physics.characteristic)(ctx->physics.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds,ctx->xmin+hx*i)); 382 PetscCall(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 PetscCall((*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 PetscCall((*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 PetscCall((*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 PetscCall(DMDAVecRestoreArray(da,Xloc,&x)); 455 PetscCall(VecRestoreArray(F,&f)); 456 PetscCall(DMDARestoreArray(da,PETSC_TRUE,&slope)); 457 PetscCall(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 PetscCall(TSGetDM(ts,&da)); 472 PetscCall(DMGetLocalVector(da,&Xloc)); 473 PetscCall(DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0)); 474 hx = (ctx->xmax-ctx->xmin)/Mx; 475 PetscCall(DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc)); 476 PetscCall(DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc)); 477 PetscCall(VecZeroEntries(F)); 478 PetscCall(DMDAVecGetArray(da,Xloc,&x)); 479 PetscCall(VecGetArray(F,&f)); 480 PetscCall(DMDAGetArray(da,PETSC_TRUE,&slope)); 481 PetscCall(DMDAGetCorners(da,&xs,0,0,&xm,0,0)); 482 PetscCall(ISGetSize(ctx->iss,&len_slow1)); 483 PetscCall(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 PetscCall((*ctx->physics.characteristic)(ctx->physics.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds,ctx->xmin+hx*i)); 498 PetscCall(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 PetscCall((*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 PetscCall((*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 PetscCall(DMDAVecRestoreArray(da,Xloc,&x)); 557 PetscCall(VecRestoreArray(F,&f)); 558 PetscCall(DMDARestoreArray(da,PETSC_TRUE,&slope)); 559 PetscCall(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 PetscCall(PetscInitialize(&argc,&argv,0,help)); 578 comm = PETSC_COMM_WORLD; 579 PetscCall(PetscMemzero(&ctx,sizeof(ctx))); 580 581 /* Register limiters to be available on the command line */ 582 PetscCall(PetscFunctionListAdd(&limiters,"upwind" ,Limit_Upwind)); 583 PetscCall(PetscFunctionListAdd(&limiters,"lax-wendroff" ,Limit_LaxWendroff)); 584 PetscCall(PetscFunctionListAdd(&limiters,"beam-warming" ,Limit_BeamWarming)); 585 PetscCall(PetscFunctionListAdd(&limiters,"fromm" ,Limit_Fromm)); 586 PetscCall(PetscFunctionListAdd(&limiters,"minmod" ,Limit_Minmod)); 587 PetscCall(PetscFunctionListAdd(&limiters,"superbee" ,Limit_Superbee)); 588 PetscCall(PetscFunctionListAdd(&limiters,"mc" ,Limit_MC)); 589 PetscCall(PetscFunctionListAdd(&limiters,"vanleer" ,Limit_VanLeer)); 590 PetscCall(PetscFunctionListAdd(&limiters,"vanalbada" ,Limit_VanAlbada)); 591 PetscCall(PetscFunctionListAdd(&limiters,"vanalbadatvd" ,Limit_VanAlbadaTVD)); 592 PetscCall(PetscFunctionListAdd(&limiters,"koren" ,Limit_Koren)); 593 PetscCall(PetscFunctionListAdd(&limiters,"korensym" ,Limit_KorenSym)); 594 PetscCall(PetscFunctionListAdd(&limiters,"koren3" ,Limit_Koren3)); 595 PetscCall(PetscFunctionListAdd(&limiters,"cada-torrilhon2" ,Limit_CadaTorrilhon2)); 596 PetscCall(PetscFunctionListAdd(&limiters,"cada-torrilhon3-r0p1",Limit_CadaTorrilhon3R0p1)); 597 PetscCall(PetscFunctionListAdd(&limiters,"cada-torrilhon3-r1" ,Limit_CadaTorrilhon3R1)); 598 PetscCall(PetscFunctionListAdd(&limiters,"cada-torrilhon3-r10" ,Limit_CadaTorrilhon3R10)); 599 PetscCall(PetscFunctionListAdd(&limiters,"cada-torrilhon3-r100",Limit_CadaTorrilhon3R100)); 600 601 /* Register physical models to be available on the command line */ 602 PetscCall(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","");PetscCall(ierr); 610 PetscCall(PetscOptionsReal("-xmin","X min","",ctx.xmin,&ctx.xmin,NULL)); 611 PetscCall(PetscOptionsReal("-xmax","X max","",ctx.xmax,&ctx.xmax,NULL)); 612 PetscCall(PetscOptionsFList("-limit","Name of flux limiter to use","",limiters,lname,lname,sizeof(lname),NULL)); 613 PetscCall(PetscOptionsInt("-draw","Draw solution vector, bitwise OR of (1=initial,2=final,4=final error)","",draw,&draw,NULL)); 614 PetscCall(PetscOptionsString("-view_final","Write final solution in ASCII MATLAB format to given file name","",final_fname,final_fname,sizeof(final_fname),&view_final)); 615 PetscCall(PetscOptionsInt("-initial","Initial condition (depends on the physics)","",ctx.initial,&ctx.initial,NULL)); 616 PetscCall(PetscOptionsBool("-simulation","Compare errors with reference solution","",ctx.simulation,&ctx.simulation,NULL)); 617 PetscCall(PetscOptionsReal("-cfl","CFL number to time step at","",ctx.cfl,&ctx.cfl,NULL)); 618 PetscCall(PetscOptionsEnum("-bc_type","Boundary condition","",FVBCTypes,(PetscEnum)ctx.bctype,(PetscEnum*)&ctx.bctype,NULL)); 619 PetscCall(PetscOptionsBool("-recursive_split","Split the domain recursively","",ctx.recursive,&ctx.recursive,NULL)); 620 ierr = PetscOptionsEnd();PetscCall(ierr); 621 622 /* Choose the limiter from the list of registered limiters */ 623 PetscCall(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 PetscCall(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 PetscCall((*r)(&ctx)); 633 } 634 635 /* Create a DMDA to manage the parallel grid */ 636 PetscCall(DMDACreate1d(comm,DM_BOUNDARY_PERIODIC,50,ctx.physics.dof,2,NULL,&da)); 637 PetscCall(DMSetFromOptions(da)); 638 PetscCall(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 PetscCall(DMDASetFieldName(da,i,ctx.physics.fieldname[i])); 643 } 644 PetscCall(DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0)); 645 PetscCall(DMDAGetCorners(da,&xs,0,0,&xm,0,0)); 646 647 /* Set coordinates of cell centers */ 648 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)); 649 650 /* Allocate work space for the Finite Volume solver (so it doesn't have to be reallocated on each function evaluation) */ 651 PetscCall(PetscMalloc4(dof*dof,&ctx.R,dof*dof,&ctx.Rinv,2*dof,&ctx.cjmpLR,1*dof,&ctx.cslope)); 652 PetscCall(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 PetscCall(DMCreateGlobalVector(da,&X)); 656 PetscCall(VecDuplicate(X,&X0)); 657 PetscCall(VecDuplicate(X,&R)); 658 659 /*-------------------------------- create index for slow parts and fast parts ----------------------------------------*/ 660 PetscCall(PetscMalloc1(xm*dof,&index_slow)); 661 PetscCall(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 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD,islow,index_slow,PETSC_COPY_VALUES,&ctx.iss)); 669 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD,ifast,index_fast,PETSC_COPY_VALUES,&ctx.isf)); 670 671 /* Create a time-stepping object */ 672 PetscCall(TSCreate(comm,&ts)); 673 PetscCall(TSSetDM(ts,da)); 674 PetscCall(TSSetRHSFunction(ts,R,FVRHSFunction,&ctx)); 675 PetscCall(TSRHSSplitSetIS(ts,"slow",ctx.iss)); 676 PetscCall(TSRHSSplitSetIS(ts,"fast",ctx.isf)); 677 PetscCall(TSRHSSplitSetRHSFunction(ts,"slow",NULL,FVRHSFunctionslow,&ctx)); 678 PetscCall(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 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD,islow,index_slow,PETSC_COPY_VALUES,&ctx.iss2)); 692 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD,ifast,index_fast,PETSC_COPY_VALUES,&ctx.isf2)); 693 694 PetscCall(TSRHSSplitGetSubTS(ts,"fast",&subts)); 695 PetscCall(TSRHSSplitSetIS(subts,"slow",ctx.iss2)); 696 PetscCall(TSRHSSplitSetIS(subts,"fast",ctx.isf2)); 697 PetscCall(TSRHSSplitSetRHSFunction(subts,"slow",NULL,FVRHSFunctionslow2,&ctx)); 698 PetscCall(TSRHSSplitSetRHSFunction(subts,"fast",NULL,FVRHSFunctionfast2,&ctx)); 699 } 700 701 /*PetscCall(TSSetType(ts,TSSSP));*/ 702 PetscCall(TSSetType(ts,TSMPRK)); 703 PetscCall(TSSetMaxTime(ts,10)); 704 PetscCall(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER)); 705 706 /* Compute initial conditions and starting time step */ 707 PetscCall(FVSample(&ctx,da,0,X0)); 708 PetscCall(FVRHSFunction(ts,0,X0,X,(void*)&ctx)); /* Initial function evaluation, only used to determine max speed */ 709 PetscCall(VecCopy(X0,X)); /* The function value was not used so we set X=X0 again */ 710 PetscCall(TSSetTimeStep(ts,ctx.cfl/ctx.cfl_idt)); 711 PetscCall(TSSetFromOptions(ts)); /* Take runtime options */ 712 PetscCall(SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD)); 713 { 714 PetscInt steps; 715 PetscScalar mass_initial,mass_final,mass_difference; 716 717 PetscCall(TSSolve(ts,X)); 718 PetscCall(TSGetSolveTime(ts,&ptime)); 719 PetscCall(TSGetStepNumber(ts,&steps)); 720 PetscCall(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 PetscCall(VecSum(X0,&mass_initial)); 725 PetscCall(VecSum(X,&mass_final)); 726 mass_difference = (ctx.xmax-ctx.xmin)/(PetscScalar)Mx*(mass_final - mass_initial); 727 PetscCall(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 PetscCall(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 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&fd)); 738 PetscCall(VecDuplicate(X0,&XR)); 739 PetscCall(VecLoad(XR,fd)); 740 PetscCall(PetscViewerDestroy(&fd)); 741 PetscCall(VecAYPX(XR,-1,X)); 742 PetscCall(VecNorm(XR,NORM_1,&nrm1)); 743 PetscCall(VecNorm(XR,NORM_INFINITY,&nrmsup)); 744 nrm1 /= Mx; 745 PetscCall(PetscPrintf(comm,"Error ||x-x_e||_1 %g ||x-x_e||_sup %g\n",(double)nrm1,(double)nrmsup)); 746 PetscCall(VecDestroy(&XR)); 747 } 748 } 749 750 PetscCall(SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD)); 751 if (draw & 0x1) PetscCall(VecView(X0,PETSC_VIEWER_DRAW_WORLD)); 752 if (draw & 0x2) PetscCall(VecView(X,PETSC_VIEWER_DRAW_WORLD)); 753 if (draw & 0x4) { 754 Vec Y; 755 PetscCall(VecDuplicate(X,&Y)); 756 PetscCall(FVSample(&ctx,da,ptime,Y)); 757 PetscCall(VecAYPX(Y,-1,X)); 758 PetscCall(VecView(Y,PETSC_VIEWER_DRAW_WORLD)); 759 PetscCall(VecDestroy(&Y)); 760 } 761 762 if (view_final) { 763 PetscViewer viewer; 764 PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD,final_fname,&viewer)); 765 PetscCall(PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB)); 766 PetscCall(VecView(X,viewer)); 767 PetscCall(PetscViewerPopFormat(viewer)); 768 PetscCall(PetscViewerDestroy(&viewer)); 769 } 770 771 /* Clean up */ 772 PetscCall((*ctx.physics.destroy)(ctx.physics.user)); 773 for (i=0; i<ctx.physics.dof; i++) PetscCall(PetscFree(ctx.physics.fieldname[i])); 774 PetscCall(PetscFree4(ctx.R,ctx.Rinv,ctx.cjmpLR,ctx.cslope)); 775 PetscCall(PetscFree3(ctx.uLR,ctx.flux,ctx.speeds)); 776 PetscCall(ISDestroy(&ctx.iss)); 777 PetscCall(ISDestroy(&ctx.isf)); 778 PetscCall(ISDestroy(&ctx.iss2)); 779 PetscCall(ISDestroy(&ctx.isf2)); 780 PetscCall(VecDestroy(&X)); 781 PetscCall(VecDestroy(&X0)); 782 PetscCall(VecDestroy(&R)); 783 PetscCall(DMDestroy(&da)); 784 PetscCall(TSDestroy(&ts)); 785 PetscCall(PetscFree(index_slow)); 786 PetscCall(PetscFree(index_fast)); 787 PetscCall(PetscFunctionListDestroy(&limiters)); 788 PetscCall(PetscFunctionListDestroy(&physics)); 789 PetscCall(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