1 /* 2 Note: 3 -hratio is the ratio between mesh size of corse grids and fine grids 4 -ts_rk_dtratio is the ratio between the large stepsize and the small stepsize 5 */ 6 7 static const char help[] = "1D periodic Finite Volume solver in slope-limiter form with semidiscrete time stepping.\n" 8 " advection - Constant coefficient scalar advection\n" 9 " u_t + (a*u)_x = 0\n" 10 " for this toy problem, we choose different meshsizes for different sub-domains (slow-fast-slow), say\n" 11 " hxs = hratio*hxf \n" 12 " where hxs and hxf are the grid spacings for coarse and fine grids respectively.\n" 13 " exact - Exact Riemann solver which usually needs to perform a Newton iteration to connect\n" 14 " the states across shocks and rarefactions\n" 15 " simulation - use reference solution which is generated by smaller time step size to be true solution,\n" 16 " also the reference solution should be generated by user and stored in a binary file.\n" 17 " characteristic - Limit the characteristic variables, this is usually preferred (default)\n" 18 "Several initial conditions can be chosen with -initial N\n\n" 19 "The problem size should be set with -da_grid_x M\n\n"; 20 21 #include <petscts.h> 22 #include <petscdm.h> 23 #include <petscdmda.h> 24 #include <petscdraw.h> 25 #include "finitevolume1d.h" 26 27 static inline PetscReal RangeMod(PetscReal a,PetscReal xmin,PetscReal xmax) { PetscReal range = xmax-xmin; return xmin +PetscFmodReal(range+PetscFmodReal(a,range),range); } 28 29 /* --------------------------------- Advection ----------------------------------- */ 30 typedef struct { 31 PetscReal a; /* advective velocity */ 32 } AdvectCtx; 33 34 static PetscErrorCode PhysicsRiemann_Advect(void *vctx,PetscInt m,const PetscScalar *uL,const PetscScalar *uR,PetscScalar *flux,PetscReal *maxspeed) 35 { 36 AdvectCtx *ctx = (AdvectCtx*)vctx; 37 PetscReal speed; 38 39 PetscFunctionBeginUser; 40 speed = ctx->a; 41 flux[0] = PetscMax(0,speed)*uL[0] + PetscMin(0,speed)*uR[0]; 42 *maxspeed = speed; 43 PetscFunctionReturn(0); 44 } 45 46 static PetscErrorCode PhysicsCharacteristic_Advect(void *vctx,PetscInt m,const PetscScalar *u,PetscScalar *X,PetscScalar *Xi,PetscReal *speeds) 47 { 48 AdvectCtx *ctx = (AdvectCtx*)vctx; 49 50 PetscFunctionBeginUser; 51 X[0] = 1.; 52 Xi[0] = 1.; 53 speeds[0] = ctx->a; 54 PetscFunctionReturn(0); 55 } 56 57 static PetscErrorCode PhysicsSample_Advect(void *vctx,PetscInt initial,FVBCType bctype,PetscReal xmin,PetscReal xmax,PetscReal t,PetscReal x,PetscReal *u) 58 { 59 AdvectCtx *ctx = (AdvectCtx*)vctx; 60 PetscReal a = ctx->a,x0; 61 62 PetscFunctionBeginUser; 63 switch (bctype) { 64 case FVBC_OUTFLOW: x0 = x-a*t; break; 65 case FVBC_PERIODIC: x0 = RangeMod(x-a*t,xmin,xmax); break; 66 default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown BCType"); 67 } 68 switch (initial) { 69 case 0: u[0] = (x0 < 0) ? 1 : -1; break; 70 case 1: u[0] = (x0 < 0) ? -1 : 1; break; 71 case 2: u[0] = (0 < x0 && x0 < 1) ? 1 : 0; break; 72 case 3: u[0] = PetscSinReal(2*PETSC_PI*x0); break; 73 case 4: u[0] = PetscAbs(x0); break; 74 case 5: u[0] = (x0 < 0 || x0 > 0.5) ? 0 : PetscSqr(PetscSinReal(2*PETSC_PI*x0)); break; 75 case 6: u[0] = (x0 < 0) ? 0 : ((x0 < 1) ? x0 : ((x0 < 2) ? 2-x0 : 0)); break; 76 case 7: u[0] = PetscPowReal(PetscSinReal(PETSC_PI*x0),10.0);break; 77 default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"unknown initial condition"); 78 } 79 PetscFunctionReturn(0); 80 } 81 82 static PetscErrorCode PhysicsCreate_Advect(FVCtx *ctx) 83 { 84 AdvectCtx *user; 85 86 PetscFunctionBeginUser; 87 PetscCall(PetscNew(&user)); 88 ctx->physics2.sample2 = PhysicsSample_Advect; 89 ctx->physics2.riemann2 = PhysicsRiemann_Advect; 90 ctx->physics2.characteristic2 = PhysicsCharacteristic_Advect; 91 ctx->physics2.destroy = PhysicsDestroy_SimpleFree; 92 ctx->physics2.user = user; 93 ctx->physics2.dof = 1; 94 PetscCall(PetscStrallocpy("u",&ctx->physics2.fieldname[0])); 95 user->a = 1; 96 PetscOptionsBegin(ctx->comm,ctx->prefix,"Options for advection",""); 97 { 98 PetscCall(PetscOptionsReal("-physics_advect_a","Speed","",user->a,&user->a,NULL)); 99 } 100 PetscOptionsEnd(); 101 PetscFunctionReturn(0); 102 } 103 104 PetscErrorCode FVSample_2WaySplit(FVCtx *ctx,DM da,PetscReal time,Vec U) 105 { 106 PetscScalar *u,*uj,xj,xi; 107 PetscInt i,j,k,dof,xs,xm,Mx; 108 const PetscInt N = 200; 109 PetscReal hs,hf; 110 111 PetscFunctionBeginUser; 112 PetscCheck(ctx->physics2.sample2,PETSC_COMM_SELF,PETSC_ERR_SUP,"Physics has not provided a sampling function"); 113 PetscCall(DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0)); 114 PetscCall(DMDAGetCorners(da,&xs,0,0,&xm,0,0)); 115 PetscCall(DMDAVecGetArray(da,U,&u)); 116 PetscCall(PetscMalloc1(dof,&uj)); 117 hs = (ctx->xmax-ctx->xmin)*3.0/8.0/ctx->sf; 118 hf = (ctx->xmax-ctx->xmin)/4.0/(ctx->fs-ctx->sf); 119 for (i=xs; i<xs+xm; i++) { 120 if (i < ctx->sf) { 121 xi = ctx->xmin+0.5*hs+i*hs; 122 /* Integrate over cell i using trapezoid rule with N points. */ 123 for (k=0; k<dof; k++) u[i*dof+k] = 0; 124 for (j=0; j<N+1; j++) { 125 xj = xi+hs*(j-N/2)/(PetscReal)N; 126 PetscCall((*ctx->physics2.sample2)(ctx->physics2.user,ctx->initial,ctx->bctype,ctx->xmin,ctx->xmax,time,xj,uj)); 127 for (k=0; k<dof; k++) u[i*dof+k] += ((j==0 || j==N) ? 0.5 : 1.0)*uj[k]/N; 128 } 129 } else if (i < ctx->fs) { 130 xi = ctx->xmin+ctx->sf*hs+0.5*hf+(i-ctx->sf)*hf; 131 /* Integrate over cell i using trapezoid rule with N points. */ 132 for (k=0; k<dof; k++) u[i*dof+k] = 0; 133 for (j=0; j<N+1; j++) { 134 xj = xi+hf*(j-N/2)/(PetscReal)N; 135 PetscCall((*ctx->physics2.sample2)(ctx->physics2.user,ctx->initial,ctx->bctype,ctx->xmin,ctx->xmax,time,xj,uj)); 136 for (k=0; k<dof; k++) u[i*dof+k] += ((j==0 || j==N) ? 0.5 : 1.0)*uj[k]/N; 137 } 138 } else { 139 xi = ctx->xmin+ctx->sf*hs+(ctx->fs-ctx->sf)*hf+0.5*hs+(i-ctx->fs)*hs; 140 /* Integrate over cell i using trapezoid rule with N points. */ 141 for (k=0; k<dof; k++) u[i*dof+k] = 0; 142 for (j=0; j<N+1; j++) { 143 xj = xi+hs*(j-N/2)/(PetscReal)N; 144 PetscCall((*ctx->physics2.sample2)(ctx->physics2.user,ctx->initial,ctx->bctype,ctx->xmin,ctx->xmax,time,xj,uj)); 145 for (k=0; k<dof; k++) u[i*dof+k] += ((j==0 || j==N) ? 0.5 : 1.0)*uj[k]/N; 146 } 147 } 148 } 149 PetscCall(DMDAVecRestoreArray(da,U,&u)); 150 PetscCall(PetscFree(uj)); 151 PetscFunctionReturn(0); 152 } 153 154 static PetscErrorCode SolutionErrorNorms_2WaySplit(FVCtx *ctx,DM da,PetscReal t,Vec X,PetscReal *nrm1) 155 { 156 Vec Y; 157 PetscInt i,Mx; 158 const PetscScalar *ptr_X,*ptr_Y; 159 PetscReal hs,hf; 160 161 PetscFunctionBeginUser; 162 PetscCall(VecGetSize(X,&Mx)); 163 PetscCall(VecDuplicate(X,&Y)); 164 PetscCall(FVSample_2WaySplit(ctx,da,t,Y)); 165 hs = (ctx->xmax-ctx->xmin)*3.0/8.0/ctx->sf; 166 hf = (ctx->xmax-ctx->xmin)/4.0/(ctx->fs-ctx->sf); 167 PetscCall(VecGetArrayRead(X,&ptr_X)); 168 PetscCall(VecGetArrayRead(Y,&ptr_Y)); 169 for (i=0; i<Mx; i++) { 170 if (i < ctx->sf || i > ctx->fs -1) *nrm1 += hs*PetscAbs(ptr_X[i]-ptr_Y[i]); 171 else *nrm1 += hf*PetscAbs(ptr_X[i]-ptr_Y[i]); 172 } 173 PetscCall(VecRestoreArrayRead(X,&ptr_X)); 174 PetscCall(VecRestoreArrayRead(Y,&ptr_Y)); 175 PetscCall(VecDestroy(&Y)); 176 PetscFunctionReturn(0); 177 } 178 179 PetscErrorCode FVRHSFunction_2WaySplit(TS ts,PetscReal time,Vec X,Vec F,void *vctx) 180 { 181 FVCtx *ctx = (FVCtx*)vctx; 182 PetscInt i,j,k,Mx,dof,xs,xm,sf = ctx->sf,fs = ctx->fs; 183 PetscReal hxf,hxs,cfl_idt = 0; 184 PetscScalar *x,*f,*slope; 185 Vec Xloc; 186 DM da; 187 188 PetscFunctionBeginUser; 189 PetscCall(TSGetDM(ts,&da)); 190 PetscCall(DMGetLocalVector(da,&Xloc)); /* Xloc contains ghost points */ 191 PetscCall(DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0)); /* Mx is the number of center points */ 192 hxs = (ctx->xmax-ctx->xmin)*3.0/8.0/ctx->sf; 193 hxf = (ctx->xmax-ctx->xmin)/4.0/(ctx->fs-ctx->sf); 194 PetscCall(DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc)); /* X is solution vector which does not contain ghost points */ 195 PetscCall(DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc)); 196 197 PetscCall(VecZeroEntries(F)); /* F is the right hand side function corresponds to center points */ 198 199 PetscCall(DMDAVecGetArray(da,Xloc,&x)); 200 PetscCall(DMDAVecGetArray(da,F,&f)); 201 PetscCall(DMDAGetArray(da,PETSC_TRUE,&slope)); /* contains ghost points */ 202 203 PetscCall(DMDAGetCorners(da,&xs,0,0,&xm,0,0)); 204 205 if (ctx->bctype == FVBC_OUTFLOW) { 206 for (i=xs-2; i<0; i++) { 207 for (j=0; j<dof; j++) x[i*dof+j] = x[j]; 208 } 209 for (i=Mx; i<xs+xm+2; i++) { 210 for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j]; 211 } 212 } 213 for (i=xs-1; i<xs+xm+1; i++) { 214 struct _LimitInfo info; 215 PetscScalar *cjmpL,*cjmpR; 216 /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */ 217 PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds)); 218 /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */ 219 PetscCall(PetscArrayzero(ctx->cjmpLR,2*dof)); 220 cjmpL = &ctx->cjmpLR[0]; 221 cjmpR = &ctx->cjmpLR[dof]; 222 for (j=0; j<dof; j++) { 223 PetscScalar jmpL,jmpR; 224 jmpL = x[(i+0)*dof+j]-x[(i-1)*dof+j]; 225 jmpR = x[(i+1)*dof+j]-x[(i+0)*dof+j]; 226 for (k=0; k<dof; k++) { 227 cjmpL[k] += ctx->Rinv[k+j*dof]*jmpL; 228 cjmpR[k] += ctx->Rinv[k+j*dof]*jmpR; 229 } 230 } 231 /* Apply limiter to the left and right characteristic jumps */ 232 info.m = dof; 233 info.hxs = hxs; 234 info.hxf = hxf; 235 (*ctx->limit2)(&info,cjmpL,cjmpR,ctx->sf,ctx->fs,i,ctx->cslope); 236 for (j=0; j<dof; j++) { 237 PetscScalar tmp = 0; 238 for (k=0; k<dof; k++) tmp += ctx->R[j+k*dof]*ctx->cslope[k]; 239 slope[i*dof+j] = tmp; 240 } 241 } 242 243 for (i=xs; i<xs+xm+1; i++) { 244 PetscReal maxspeed; 245 PetscScalar *uL,*uR; 246 uL = &ctx->uLR[0]; 247 uR = &ctx->uLR[dof]; 248 if (i < sf) { /* slow region */ 249 for (j=0; j<dof; j++) { 250 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2; 251 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2; 252 } 253 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed)); 254 if (i > xs) { 255 for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hxs; 256 } 257 if (i < xs+xm) { 258 for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hxs; 259 } 260 } else if (i == sf) { /* interface between the slow region and the fast region */ 261 for (j=0; j<dof; j++) { 262 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2; 263 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxf/2; 264 } 265 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed)); 266 if (i > xs) { 267 for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hxs; 268 } 269 if (i < xs+xm) { 270 for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hxf; 271 } 272 } else if (i > sf && i < fs) { /* fast region */ 273 for (j=0; j<dof; j++) { 274 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxf/2; 275 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxf/2; 276 } 277 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed)); 278 if (i > xs) { 279 for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hxf; 280 } 281 if (i < xs+xm) { 282 for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hxf; 283 } 284 } else if (i == fs) { /* interface between the fast region and the slow region */ 285 for (j=0; j<dof; j++) { 286 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxf/2; 287 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2; 288 } 289 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed)); 290 if (i > xs) { 291 for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hxf; 292 } 293 if (i < xs+xm) { 294 for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hxs; 295 } 296 } else { /* slow region */ 297 for (j=0; j<dof; j++) { 298 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2; 299 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2; 300 } 301 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed)); 302 cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hxs)); /* Max allowable value of 1/Delta t */ 303 if (i > xs) { 304 for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hxs; 305 } 306 if (i < xs+xm) { 307 for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hxs; 308 } 309 } 310 } 311 PetscCall(DMDAVecRestoreArray(da,Xloc,&x)); 312 PetscCall(DMDAVecRestoreArray(da,F,&f)); 313 PetscCall(DMDARestoreArray(da,PETSC_TRUE,&slope)); 314 PetscCall(DMRestoreLocalVector(da,&Xloc)); 315 PetscCallMPI(MPI_Allreduce(&cfl_idt,&ctx->cfl_idt,1,MPIU_SCALAR,MPIU_MAX,PetscObjectComm((PetscObject)da))); 316 if (0) { 317 /* We need to a way to inform the TS of a CFL constraint, this is a debugging fragment */ 318 PetscReal dt,tnow; 319 PetscCall(TSGetTimeStep(ts,&dt)); 320 PetscCall(TSGetTime(ts,&tnow)); 321 if (dt > 0.5/ctx->cfl_idt) { 322 if (1) { 323 PetscCall(PetscPrintf(ctx->comm,"Stability constraint exceeded at t=%g, dt %g > %g\n",(double)tnow,(double)dt,(double)(0.5/ctx->cfl_idt))); 324 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Stability constraint exceeded, %g > %g",(double)dt,(double)(ctx->cfl/ctx->cfl_idt)); 325 } 326 } 327 PetscFunctionReturn(0); 328 } 329 330 /* --------------------------------- Finite Volume Solver for slow components ----------------------------------- */ 331 PetscErrorCode FVRHSFunctionslow_2WaySplit(TS ts,PetscReal time,Vec X,Vec F,void *vctx) 332 { 333 FVCtx *ctx = (FVCtx*)vctx; 334 PetscInt i,j,k,Mx,dof,xs,xm,islow = 0,sf = ctx->sf,fs = ctx->fs,lsbwidth = ctx->lsbwidth,rsbwidth = ctx->rsbwidth; 335 PetscReal hxs,hxf,cfl_idt = 0; 336 PetscScalar *x,*f,*slope; 337 Vec Xloc; 338 DM da; 339 340 PetscFunctionBeginUser; 341 PetscCall(TSGetDM(ts,&da)); 342 PetscCall(DMGetLocalVector(da,&Xloc)); 343 PetscCall(DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0)); 344 hxs = (ctx->xmax-ctx->xmin)*3.0/8.0/ctx->sf; 345 hxf = (ctx->xmax-ctx->xmin)/4.0/(ctx->fs-ctx->sf); 346 PetscCall(DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc)); 347 PetscCall(DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc)); 348 PetscCall(VecZeroEntries(F)); 349 PetscCall(DMDAVecGetArray(da,Xloc,&x)); 350 PetscCall(VecGetArray(F,&f)); 351 PetscCall(DMDAGetArray(da,PETSC_TRUE,&slope)); 352 PetscCall(DMDAGetCorners(da,&xs,0,0,&xm,0,0)); 353 354 if (ctx->bctype == FVBC_OUTFLOW) { 355 for (i=xs-2; i<0; i++) { 356 for (j=0; j<dof; j++) x[i*dof+j] = x[j]; 357 } 358 for (i=Mx; i<xs+xm+2; i++) { 359 for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j]; 360 } 361 } 362 for (i=xs-1; i<xs+xm+1; i++) { 363 struct _LimitInfo info; 364 PetscScalar *cjmpL,*cjmpR; 365 if (i < sf-lsbwidth+1 || i > fs+rsbwidth-2) { /* slow components and the first and last fast components */ 366 /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */ 367 PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds)); 368 /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */ 369 PetscCall(PetscArrayzero(ctx->cjmpLR,2*dof)); 370 cjmpL = &ctx->cjmpLR[0]; 371 cjmpR = &ctx->cjmpLR[dof]; 372 for (j=0; j<dof; j++) { 373 PetscScalar jmpL,jmpR; 374 jmpL = x[(i+0)*dof+j]-x[(i-1)*dof+j]; 375 jmpR = x[(i+1)*dof+j]-x[(i+0)*dof+j]; 376 for (k=0; k<dof; k++) { 377 cjmpL[k] += ctx->Rinv[k+j*dof]*jmpL; 378 cjmpR[k] += ctx->Rinv[k+j*dof]*jmpR; 379 } 380 } 381 /* Apply limiter to the left and right characteristic jumps */ 382 info.m = dof; 383 info.hxs = hxs; 384 info.hxf = hxf; 385 (*ctx->limit2)(&info,cjmpL,cjmpR,ctx->sf,ctx->fs,i,ctx->cslope); 386 for (j=0; j<dof; j++) { 387 PetscScalar tmp = 0; 388 for (k=0; k<dof; k++) tmp += ctx->R[j+k*dof]*ctx->cslope[k]; 389 slope[i*dof+j] = tmp; 390 } 391 } 392 } 393 394 for (i=xs; i<xs+xm+1; i++) { 395 PetscReal maxspeed; 396 PetscScalar *uL,*uR; 397 uL = &ctx->uLR[0]; 398 uR = &ctx->uLR[dof]; 399 if (i < sf-lsbwidth) { /* slow region */ 400 for (j=0; j<dof; j++) { 401 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2; 402 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2; 403 } 404 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed)); 405 cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hxs)); /* Max allowable value of 1/Delta t */ 406 if (i > xs) { 407 for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hxs; 408 } 409 if (i < xs+xm) { 410 for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hxs; 411 islow++; 412 } 413 } 414 if (i == sf-lsbwidth) { /* interface between the slow region and the fast region */ 415 for (j=0; j<dof; j++) { 416 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2; 417 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2; 418 } 419 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed)); 420 if (i > xs) { 421 for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hxs; 422 } 423 } 424 if (i == fs+rsbwidth) { /* slow region */ 425 for (j=0; j<dof; j++) { 426 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2; 427 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2; 428 } 429 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed)); 430 if (i < xs+xm) { 431 for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hxs; 432 islow++; 433 } 434 } 435 if (i > fs+rsbwidth) { /* slow region */ 436 for (j=0; j<dof; j++) { 437 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2; 438 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2; 439 } 440 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed)); 441 if (i > xs) { 442 for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hxs; 443 } 444 if (i < xs+xm) { 445 for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hxs; 446 islow++; 447 } 448 } 449 } 450 PetscCall(DMDAVecRestoreArray(da,Xloc,&x)); 451 PetscCall(VecRestoreArray(F,&f)); 452 PetscCall(DMDARestoreArray(da,PETSC_TRUE,&slope)); 453 PetscCall(DMRestoreLocalVector(da,&Xloc)); 454 PetscCallMPI(MPI_Allreduce(&cfl_idt,&ctx->cfl_idt,1,MPIU_SCALAR,MPIU_MAX,PetscObjectComm((PetscObject)da))); 455 PetscFunctionReturn(0); 456 } 457 458 PetscErrorCode FVRHSFunctionslowbuffer_2WaySplit(TS ts,PetscReal time,Vec X,Vec F,void *vctx) 459 { 460 FVCtx *ctx = (FVCtx*)vctx; 461 PetscInt i,j,k,Mx,dof,xs,xm,islow = 0,sf = ctx->sf,fs = ctx->fs,lsbwidth = ctx->lsbwidth,rsbwidth = ctx->rsbwidth; 462 PetscReal hxs,hxf; 463 PetscScalar *x,*f,*slope; 464 Vec Xloc; 465 DM da; 466 467 PetscFunctionBeginUser; 468 PetscCall(TSGetDM(ts,&da)); 469 PetscCall(DMGetLocalVector(da,&Xloc)); 470 PetscCall(DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0)); 471 hxs = (ctx->xmax-ctx->xmin)*3.0/8.0/ctx->sf; 472 hxf = (ctx->xmax-ctx->xmin)/4.0/(ctx->fs-ctx->sf); 473 PetscCall(DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc)); 474 PetscCall(DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc)); 475 PetscCall(VecZeroEntries(F)); 476 PetscCall(DMDAVecGetArray(da,Xloc,&x)); 477 PetscCall(VecGetArray(F,&f)); 478 PetscCall(DMDAGetArray(da,PETSC_TRUE,&slope)); 479 PetscCall(DMDAGetCorners(da,&xs,0,0,&xm,0,0)); 480 481 if (ctx->bctype == FVBC_OUTFLOW) { 482 for (i=xs-2; i<0; i++) { 483 for (j=0; j<dof; j++) x[i*dof+j] = x[j]; 484 } 485 for (i=Mx; i<xs+xm+2; i++) { 486 for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j]; 487 } 488 } 489 for (i=xs-1; i<xs+xm+1; i++) { 490 struct _LimitInfo info; 491 PetscScalar *cjmpL,*cjmpR; 492 if ((i > sf-lsbwidth-2 && i < sf+1) || (i > fs-2 && i < fs+rsbwidth+1)) { 493 /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */ 494 PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds)); 495 /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */ 496 PetscCall(PetscArrayzero(ctx->cjmpLR,2*dof)); 497 cjmpL = &ctx->cjmpLR[0]; 498 cjmpR = &ctx->cjmpLR[dof]; 499 for (j=0; j<dof; j++) { 500 PetscScalar jmpL,jmpR; 501 jmpL = x[(i+0)*dof+j]-x[(i-1)*dof+j]; 502 jmpR = x[(i+1)*dof+j]-x[(i+0)*dof+j]; 503 for (k=0; k<dof; k++) { 504 cjmpL[k] += ctx->Rinv[k+j*dof]*jmpL; 505 cjmpR[k] += ctx->Rinv[k+j*dof]*jmpR; 506 } 507 } 508 /* Apply limiter to the left and right characteristic jumps */ 509 info.m = dof; 510 info.hxs = hxs; 511 info.hxf = hxf; 512 (*ctx->limit2)(&info,cjmpL,cjmpR,ctx->sf,ctx->fs,i,ctx->cslope); 513 for (j=0; j<dof; j++) { 514 PetscScalar tmp = 0; 515 for (k=0; k<dof; k++) tmp += ctx->R[j+k*dof]*ctx->cslope[k]; 516 slope[i*dof+j] = tmp; 517 } 518 } 519 } 520 521 for (i=xs; i<xs+xm+1; i++) { 522 PetscReal maxspeed; 523 PetscScalar *uL,*uR; 524 uL = &ctx->uLR[0]; 525 uR = &ctx->uLR[dof]; 526 if (i == sf-lsbwidth) { 527 for (j=0; j<dof; j++) { 528 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2; 529 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2; 530 } 531 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed)); 532 if (i < xs+xm) { 533 for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hxs; 534 islow++; 535 } 536 } 537 if (i > sf-lsbwidth && i < sf) { 538 for (j=0; j<dof; j++) { 539 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2; 540 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2; 541 } 542 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed)); 543 if (i > xs) { 544 for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hxs; 545 } 546 if (i < xs+xm) { 547 for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hxs; 548 islow++; 549 } 550 } 551 if (i == sf) { /* interface between the slow region and the fast region */ 552 for (j=0; j<dof; j++) { 553 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2; 554 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxf/2; 555 } 556 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed)); 557 if (i > xs) { 558 for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hxs; 559 } 560 } 561 if (i == fs) { /* interface between the fast region and the slow region */ 562 for (j=0; j<dof; j++) { 563 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxf/2; 564 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2; 565 } 566 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed)); 567 if (i < xs+xm) { 568 for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hxs; 569 islow++; 570 } 571 } 572 if (i > fs && i < fs+rsbwidth) { 573 for (j=0; j<dof; j++) { 574 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2; 575 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2; 576 } 577 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed)); 578 if (i > xs) { 579 for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hxs; 580 } 581 if (i < xs+xm) { 582 for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hxs; 583 islow++; 584 } 585 } 586 if (i == fs+rsbwidth) { 587 for (j=0; j<dof; j++) { 588 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2; 589 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2; 590 } 591 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed)); 592 if (i > xs) { 593 for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hxs; 594 } 595 } 596 } 597 PetscCall(DMDAVecRestoreArray(da,Xloc,&x)); 598 PetscCall(VecRestoreArray(F,&f)); 599 PetscCall(DMDARestoreArray(da,PETSC_TRUE,&slope)); 600 PetscCall(DMRestoreLocalVector(da,&Xloc)); 601 PetscFunctionReturn(0); 602 } 603 604 /* --------------------------------- Finite Volume Solver for fast parts ----------------------------------- */ 605 PetscErrorCode FVRHSFunctionfast_2WaySplit(TS ts,PetscReal time,Vec X,Vec F,void *vctx) 606 { 607 FVCtx *ctx = (FVCtx*)vctx; 608 PetscInt i,j,k,Mx,dof,xs,xm,ifast = 0,sf = ctx->sf,fs = ctx->fs; 609 PetscReal hxs,hxf; 610 PetscScalar *x,*f,*slope; 611 Vec Xloc; 612 DM da; 613 614 PetscFunctionBeginUser; 615 PetscCall(TSGetDM(ts,&da)); 616 PetscCall(DMGetLocalVector(da,&Xloc)); 617 PetscCall(DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0)); 618 hxs = (ctx->xmax-ctx->xmin)*3.0/8.0/ctx->sf; 619 hxf = (ctx->xmax-ctx->xmin)/4.0/(ctx->fs-ctx->sf); 620 PetscCall(DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc)); 621 PetscCall(DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc)); 622 PetscCall(VecZeroEntries(F)); 623 PetscCall(DMDAVecGetArray(da,Xloc,&x)); 624 PetscCall(VecGetArray(F,&f)); 625 PetscCall(DMDAGetArray(da,PETSC_TRUE,&slope)); 626 PetscCall(DMDAGetCorners(da,&xs,0,0,&xm,0,0)); 627 628 if (ctx->bctype == FVBC_OUTFLOW) { 629 for (i=xs-2; i<0; i++) { 630 for (j=0; j<dof; j++) x[i*dof+j] = x[j]; 631 } 632 for (i=Mx; i<xs+xm+2; i++) { 633 for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j]; 634 } 635 } 636 for (i=xs-1; i<xs+xm+1; i++) { 637 struct _LimitInfo info; 638 PetscScalar *cjmpL,*cjmpR; 639 if (i > sf-2 && i < fs+1) { 640 PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds)); 641 PetscCall(PetscArrayzero(ctx->cjmpLR,2*dof)); 642 cjmpL = &ctx->cjmpLR[0]; 643 cjmpR = &ctx->cjmpLR[dof]; 644 for (j=0; j<dof; j++) { 645 PetscScalar jmpL,jmpR; 646 jmpL = x[(i+0)*dof+j]-x[(i-1)*dof+j]; 647 jmpR = x[(i+1)*dof+j]-x[(i+0)*dof+j]; 648 for (k=0; k<dof; k++) { 649 cjmpL[k] += ctx->Rinv[k+j*dof]*jmpL; 650 cjmpR[k] += ctx->Rinv[k+j*dof]*jmpR; 651 } 652 } 653 /* Apply limiter to the left and right characteristic jumps */ 654 info.m = dof; 655 info.hxs = hxs; 656 info.hxf = hxf; 657 (*ctx->limit2)(&info,cjmpL,cjmpR,ctx->sf,ctx->fs,i,ctx->cslope); 658 for (j=0; j<dof; j++) { 659 PetscScalar tmp = 0; 660 for (k=0; k<dof; k++) tmp += ctx->R[j+k*dof]*ctx->cslope[k]; 661 slope[i*dof+j] = tmp; 662 } 663 } 664 } 665 666 for (i=xs; i<xs+xm+1; i++) { 667 PetscReal maxspeed; 668 PetscScalar *uL,*uR; 669 uL = &ctx->uLR[0]; 670 uR = &ctx->uLR[dof]; 671 if (i == sf) { /* interface between the slow region and the fast region */ 672 for (j=0; j<dof; j++) { 673 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxs/2; 674 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxf/2; 675 } 676 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed)); 677 if (i < xs+xm) { 678 for (j=0; j<dof; j++) f[ifast*dof+j] += ctx->flux[j]/hxf; 679 ifast++; 680 } 681 } 682 if (i > sf && i < fs) { /* fast region */ 683 for (j=0; j<dof; j++) { 684 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxf/2; 685 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxf/2; 686 } 687 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed)); 688 if (i > xs) { 689 for (j=0; j<dof; j++) f[(ifast-1)*dof+j] -= ctx->flux[j]/hxf; 690 } 691 if (i < xs+xm) { 692 for (j=0; j<dof; j++) f[ifast*dof+j] += ctx->flux[j]/hxf; 693 ifast++; 694 } 695 } 696 if (i == fs) { /* interface between the fast region and the slow region */ 697 for (j=0; j<dof; j++) { 698 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hxf/2; 699 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hxs/2; 700 } 701 PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user,dof,uL,uR,ctx->flux,&maxspeed)); 702 if (i > xs) { 703 for (j=0; j<dof; j++) f[(ifast-1)*dof+j] -= ctx->flux[j]/hxf; 704 } 705 } 706 } 707 PetscCall(DMDAVecRestoreArray(da,Xloc,&x)); 708 PetscCall(VecRestoreArray(F,&f)); 709 PetscCall(DMDARestoreArray(da,PETSC_TRUE,&slope)); 710 PetscCall(DMRestoreLocalVector(da,&Xloc)); 711 PetscFunctionReturn(0); 712 } 713 714 int main(int argc,char *argv[]) 715 { 716 char lname[256] = "mc",physname[256] = "advect",final_fname[256] = "solution.m"; 717 PetscFunctionList limiters = 0,physics = 0; 718 MPI_Comm comm; 719 TS ts; 720 DM da; 721 Vec X,X0,R; 722 FVCtx ctx; 723 PetscInt i,k,dof,xs,xm,Mx,draw = 0,count_slow,count_fast,islow = 0,ifast =0,islowbuffer = 0,*index_slow,*index_fast,*index_slowbuffer; 724 PetscBool view_final = PETSC_FALSE; 725 PetscReal ptime; 726 727 PetscCall(PetscInitialize(&argc,&argv,0,help)); 728 comm = PETSC_COMM_WORLD; 729 PetscCall(PetscMemzero(&ctx,sizeof(ctx))); 730 731 /* Register limiters to be available on the command line */ 732 PetscCall(PetscFunctionListAdd(&limiters,"upwind" ,Limit2_Upwind)); 733 PetscCall(PetscFunctionListAdd(&limiters,"lax-wendroff" ,Limit2_LaxWendroff)); 734 PetscCall(PetscFunctionListAdd(&limiters,"beam-warming" ,Limit2_BeamWarming)); 735 PetscCall(PetscFunctionListAdd(&limiters,"fromm" ,Limit2_Fromm)); 736 PetscCall(PetscFunctionListAdd(&limiters,"minmod" ,Limit2_Minmod)); 737 PetscCall(PetscFunctionListAdd(&limiters,"superbee" ,Limit2_Superbee)); 738 PetscCall(PetscFunctionListAdd(&limiters,"mc" ,Limit2_MC)); 739 PetscCall(PetscFunctionListAdd(&limiters,"koren3" ,Limit2_Koren3)); 740 741 /* Register physical models to be available on the command line */ 742 PetscCall(PetscFunctionListAdd(&physics,"advect" ,PhysicsCreate_Advect)); 743 744 ctx.comm = comm; 745 ctx.cfl = 0.9; 746 ctx.bctype = FVBC_PERIODIC; 747 ctx.xmin = -1.0; 748 ctx.xmax = 1.0; 749 PetscOptionsBegin(comm,NULL,"Finite Volume solver options",""); 750 PetscCall(PetscOptionsReal("-xmin","X min","",ctx.xmin,&ctx.xmin,NULL)); 751 PetscCall(PetscOptionsReal("-xmax","X max","",ctx.xmax,&ctx.xmax,NULL)); 752 PetscCall(PetscOptionsFList("-limit","Name of flux imiter to use","",limiters,lname,lname,sizeof(lname),NULL)); 753 PetscCall(PetscOptionsInt("-draw","Draw solution vector, bitwise OR of (1=initial,2=final,4=final error)","",draw,&draw,NULL)); 754 PetscCall(PetscOptionsString("-view_final","Write final solution in ASCII MATLAB format to given file name","",final_fname,final_fname,sizeof(final_fname),&view_final)); 755 PetscCall(PetscOptionsInt("-initial","Initial condition (depends on the physics)","",ctx.initial,&ctx.initial,NULL)); 756 PetscCall(PetscOptionsBool("-exact","Compare errors with exact solution","",ctx.exact,&ctx.exact,NULL)); 757 PetscCall(PetscOptionsBool("-simulation","Compare errors with reference solution","",ctx.simulation,&ctx.simulation,NULL)); 758 PetscCall(PetscOptionsReal("-cfl","CFL number to time step at","",ctx.cfl,&ctx.cfl,NULL)); 759 PetscCall(PetscOptionsEnum("-bc_type","Boundary condition","",FVBCTypes,(PetscEnum)ctx.bctype,(PetscEnum*)&ctx.bctype,NULL)); 760 PetscCall(PetscOptionsInt("-hratio","Spacing ratio","",ctx.hratio,&ctx.hratio,NULL)); 761 PetscOptionsEnd(); 762 763 /* Choose the limiter from the list of registered limiters */ 764 PetscCall(PetscFunctionListFind(limiters,lname,&ctx.limit2)); 765 PetscCheck(ctx.limit2,PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Limiter '%s' not found",lname); 766 767 /* Choose the physics from the list of registered models */ 768 { 769 PetscErrorCode (*r)(FVCtx*); 770 PetscCall(PetscFunctionListFind(physics,physname,&r)); 771 PetscCheck(r,PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Physics '%s' not found",physname); 772 /* Create the physics, will set the number of fields and their names */ 773 PetscCall((*r)(&ctx)); 774 } 775 776 /* Create a DMDA to manage the parallel grid */ 777 PetscCall(DMDACreate1d(comm,DM_BOUNDARY_PERIODIC,50,ctx.physics2.dof,2,NULL,&da)); 778 PetscCall(DMSetFromOptions(da)); 779 PetscCall(DMSetUp(da)); 780 /* Inform the DMDA of the field names provided by the physics. */ 781 /* The names will be shown in the title bars when run with -ts_monitor_draw_solution */ 782 for (i=0; i<ctx.physics2.dof; i++) { 783 PetscCall(DMDASetFieldName(da,i,ctx.physics2.fieldname[i])); 784 } 785 PetscCall(DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0)); 786 PetscCall(DMDAGetCorners(da,&xs,0,0,&xm,0,0)); 787 788 /* Set coordinates of cell centers */ 789 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)); 790 791 /* Allocate work space for the Finite Volume solver (so it doesn't have to be reallocated on each function evaluation) */ 792 PetscCall(PetscMalloc4(dof*dof,&ctx.R,dof*dof,&ctx.Rinv,2*dof,&ctx.cjmpLR,1*dof,&ctx.cslope)); 793 PetscCall(PetscMalloc3(2*dof,&ctx.uLR,dof,&ctx.flux,dof,&ctx.speeds)); 794 795 /* Create a vector to store the solution and to save the initial state */ 796 PetscCall(DMCreateGlobalVector(da,&X)); 797 PetscCall(VecDuplicate(X,&X0)); 798 PetscCall(VecDuplicate(X,&R)); 799 800 /* create index for slow parts and fast parts, 801 count_slow + count_fast = Mx, counts_slow*hs = 0.5, counts_fast*hf = 0.5 */ 802 count_slow = Mx/(1.0+ctx.hratio/3.0); 803 PetscCheck(count_slow%2 == 0,PETSC_COMM_WORLD,PETSC_ERR_USER,"Please adjust grid size Mx (-da_grid_x) and hratio (-hratio) so that Mx/(1+hartio/3) is even"); 804 count_fast = Mx-count_slow; 805 ctx.sf = count_slow/2; 806 ctx.fs = ctx.sf+count_fast; 807 PetscCall(PetscMalloc1(xm*dof,&index_slow)); 808 PetscCall(PetscMalloc1(xm*dof,&index_fast)); 809 PetscCall(PetscMalloc1(6*dof,&index_slowbuffer)); 810 if (((AdvectCtx*)ctx.physics2.user)->a > 0) { 811 ctx.lsbwidth = 2; 812 ctx.rsbwidth = 4; 813 } else { 814 ctx.lsbwidth = 4; 815 ctx.rsbwidth = 2; 816 } 817 for (i=xs; i<xs+xm; i++) { 818 if (i < ctx.sf-ctx.lsbwidth || i > ctx.fs+ctx.rsbwidth-1) 819 for (k=0; k<dof; k++) index_slow[islow++] = i*dof+k; 820 else if ((i >= ctx.sf-ctx.lsbwidth && i < ctx.sf) || (i > ctx.fs-1 && i <= ctx.fs+ctx.rsbwidth-1)) 821 for (k=0; k<dof; k++) index_slowbuffer[islowbuffer++] = i*dof+k; 822 else 823 for (k=0; k<dof; k++) index_fast[ifast++] = i*dof+k; 824 } 825 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD,islow,index_slow,PETSC_COPY_VALUES,&ctx.iss)); 826 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD,ifast,index_fast,PETSC_COPY_VALUES,&ctx.isf)); 827 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD,islowbuffer,index_slowbuffer,PETSC_COPY_VALUES,&ctx.issb)); 828 829 /* Create a time-stepping object */ 830 PetscCall(TSCreate(comm,&ts)); 831 PetscCall(TSSetDM(ts,da)); 832 PetscCall(TSSetRHSFunction(ts,R,FVRHSFunction_2WaySplit,&ctx)); 833 PetscCall(TSRHSSplitSetIS(ts,"slow",ctx.iss)); 834 PetscCall(TSRHSSplitSetIS(ts,"slowbuffer",ctx.issb)); 835 PetscCall(TSRHSSplitSetIS(ts,"fast",ctx.isf)); 836 PetscCall(TSRHSSplitSetRHSFunction(ts,"slow",NULL,FVRHSFunctionslow_2WaySplit,&ctx)); 837 PetscCall(TSRHSSplitSetRHSFunction(ts,"fast",NULL,FVRHSFunctionfast_2WaySplit,&ctx)); 838 PetscCall(TSRHSSplitSetRHSFunction(ts,"slowbuffer",NULL,FVRHSFunctionslowbuffer_2WaySplit,&ctx)); 839 840 PetscCall(TSSetType(ts,TSSSP)); 841 /*PetscCall(TSSetType(ts,TSMPRK));*/ 842 PetscCall(TSSetMaxTime(ts,10)); 843 PetscCall(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER)); 844 845 /* Compute initial conditions and starting time step */ 846 PetscCall(FVSample_2WaySplit(&ctx,da,0,X0)); 847 PetscCall(FVRHSFunction_2WaySplit(ts,0,X0,X,(void*)&ctx)); /* Initial function evaluation, only used to determine max speed */ 848 PetscCall(VecCopy(X0,X)); /* The function value was not used so we set X=X0 again */ 849 PetscCall(TSSetTimeStep(ts,ctx.cfl/ctx.cfl_idt)); 850 PetscCall(TSSetFromOptions(ts)); /* Take runtime options */ 851 PetscCall(SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD)); 852 { 853 PetscInt steps; 854 PetscScalar mass_initial,mass_final,mass_difference,mass_differenceg; 855 const PetscScalar *ptr_X,*ptr_X0; 856 const PetscReal hs = (ctx.xmax-ctx.xmin)*3.0/4.0/count_slow; 857 const PetscReal hf = (ctx.xmax-ctx.xmin)/4.0/count_fast; 858 859 PetscCall(TSSolve(ts,X)); 860 PetscCall(TSGetSolveTime(ts,&ptime)); 861 PetscCall(TSGetStepNumber(ts,&steps)); 862 /* calculate the total mass at initial time and final time */ 863 mass_initial = 0.0; 864 mass_final = 0.0; 865 PetscCall(DMDAVecGetArrayRead(da,X0,(void*)&ptr_X0)); 866 PetscCall(DMDAVecGetArrayRead(da,X,(void*)&ptr_X)); 867 for (i=xs;i<xs+xm;i++) { 868 if (i < ctx.sf || i > ctx.fs-1) { 869 for (k=0; k<dof; k++) { 870 mass_initial = mass_initial + hs*ptr_X0[i*dof+k]; 871 mass_final = mass_final + hs*ptr_X[i*dof+k]; 872 } 873 } else { 874 for (k=0; k<dof; k++) { 875 mass_initial = mass_initial + hf*ptr_X0[i*dof+k]; 876 mass_final = mass_final + hf*ptr_X[i*dof+k]; 877 } 878 } 879 } 880 PetscCall(DMDAVecRestoreArrayRead(da,X0,(void*)&ptr_X0)); 881 PetscCall(DMDAVecRestoreArrayRead(da,X,(void*)&ptr_X)); 882 mass_difference = mass_final - mass_initial; 883 PetscCallMPI(MPI_Allreduce(&mass_difference,&mass_differenceg,1,MPIU_SCALAR,MPIU_SUM,comm)); 884 PetscCall(PetscPrintf(comm,"Mass difference %g\n",(double)mass_differenceg)); 885 PetscCall(PetscPrintf(comm,"Final time %g, steps %" PetscInt_FMT "\n",(double)ptime,steps)); 886 PetscCall(PetscPrintf(comm,"Maximum allowable stepsize according to CFL %g\n",(double)(1.0/ctx.cfl_idt))); 887 if (ctx.exact) { 888 PetscReal nrm1=0; 889 PetscCall(SolutionErrorNorms_2WaySplit(&ctx,da,ptime,X,&nrm1)); 890 PetscCall(PetscPrintf(comm,"Error ||x-x_e||_1 %g\n",(double)nrm1)); 891 } 892 if (ctx.simulation) { 893 PetscReal nrm1=0; 894 PetscViewer fd; 895 char filename[PETSC_MAX_PATH_LEN] = "binaryoutput"; 896 Vec XR; 897 PetscBool flg; 898 const PetscScalar *ptr_XR; 899 PetscCall(PetscOptionsGetString(NULL,NULL,"-f",filename,sizeof(filename),&flg)); 900 PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_USER,"Must indicate binary file with the -f option"); 901 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&fd)); 902 PetscCall(VecDuplicate(X0,&XR)); 903 PetscCall(VecLoad(XR,fd)); 904 PetscCall(PetscViewerDestroy(&fd)); 905 PetscCall(VecGetArrayRead(X,&ptr_X)); 906 PetscCall(VecGetArrayRead(XR,&ptr_XR)); 907 for (i=xs;i<xs+xm;i++) { 908 if (i < ctx.sf || i > ctx.fs-1) 909 for (k=0; k<dof; k++) nrm1 = nrm1 + hs*PetscAbs(ptr_X[i*dof+k]-ptr_XR[i*dof+k]); 910 else 911 for (k=0; k<dof; k++) nrm1 = nrm1 + hf*PetscAbs(ptr_X[i*dof+k]-ptr_XR[i*dof+k]); 912 } 913 PetscCall(VecRestoreArrayRead(X,&ptr_X)); 914 PetscCall(VecRestoreArrayRead(XR,&ptr_XR)); 915 PetscCall(PetscPrintf(comm,"Error ||x-x_e||_1 %g\n",(double)nrm1)); 916 PetscCall(VecDestroy(&XR)); 917 } 918 } 919 920 PetscCall(SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD)); 921 if (draw & 0x1) PetscCall(VecView(X0,PETSC_VIEWER_DRAW_WORLD)); 922 if (draw & 0x2) PetscCall(VecView(X,PETSC_VIEWER_DRAW_WORLD)); 923 if (draw & 0x4) { 924 Vec Y; 925 PetscCall(VecDuplicate(X,&Y)); 926 PetscCall(FVSample_2WaySplit(&ctx,da,ptime,Y)); 927 PetscCall(VecAYPX(Y,-1,X)); 928 PetscCall(VecView(Y,PETSC_VIEWER_DRAW_WORLD)); 929 PetscCall(VecDestroy(&Y)); 930 } 931 932 if (view_final) { 933 PetscViewer viewer; 934 PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD,final_fname,&viewer)); 935 PetscCall(PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB)); 936 PetscCall(VecView(X,viewer)); 937 PetscCall(PetscViewerPopFormat(viewer)); 938 PetscCall(PetscViewerDestroy(&viewer)); 939 } 940 941 /* Clean up */ 942 PetscCall((*ctx.physics2.destroy)(ctx.physics2.user)); 943 for (i=0; i<ctx.physics2.dof; i++) PetscCall(PetscFree(ctx.physics2.fieldname[i])); 944 PetscCall(PetscFree4(ctx.R,ctx.Rinv,ctx.cjmpLR,ctx.cslope)); 945 PetscCall(PetscFree3(ctx.uLR,ctx.flux,ctx.speeds)); 946 PetscCall(VecDestroy(&X)); 947 PetscCall(VecDestroy(&X0)); 948 PetscCall(VecDestroy(&R)); 949 PetscCall(DMDestroy(&da)); 950 PetscCall(TSDestroy(&ts)); 951 PetscCall(ISDestroy(&ctx.iss)); 952 PetscCall(ISDestroy(&ctx.isf)); 953 PetscCall(ISDestroy(&ctx.issb)); 954 PetscCall(PetscFree(index_slow)); 955 PetscCall(PetscFree(index_fast)); 956 PetscCall(PetscFree(index_slowbuffer)); 957 PetscCall(PetscFunctionListDestroy(&limiters)); 958 PetscCall(PetscFunctionListDestroy(&physics)); 959 PetscCall(PetscFinalize()); 960 return 0; 961 } 962 963 /*TEST 964 965 build: 966 requires: !complex 967 depends: finitevolume1d.c 968 969 test: 970 suffix: 1 971 args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -limit mc -ts_dt 0.025 -ts_max_steps 24 -ts_type mprk -ts_mprk_type 2a22 972 973 test: 974 suffix: 2 975 args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -limit mc -ts_dt 0.025 -ts_max_steps 24 -ts_type mprk -ts_mprk_type 2a22 -ts_use_splitrhsfunction 0 976 output_file: output/ex6_1.out 977 978 TEST*/ 979