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