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