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));CHKERRMPI(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 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); 289 } 290 } 291 ierr = PetscFree4(r,min,alpha,gamma);CHKERRQ(ierr); 292 PetscFunctionReturn(0); 293 } 294 295 static PetscErrorCode FVRHSFunctionslow(TS ts,PetscReal time,Vec X,Vec F,void *vctx) 296 { 297 FVCtx *ctx = (FVCtx*)vctx; 298 PetscErrorCode ierr; 299 PetscInt i,j,Mx,dof,xs,xm,islow = 0,sf = ctx->sf,fs = ctx->fs; 300 PetscReal hf,hs; 301 PetscScalar *x,*f,*r,*min,*alpha,*gamma; 302 Vec Xloc; 303 DM da; 304 305 PetscFunctionBeginUser; 306 ierr = TSGetDM(ts,&da);CHKERRQ(ierr); 307 ierr = DMGetLocalVector(da,&Xloc);CHKERRQ(ierr); /* Xloc contains ghost points */ 308 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 */ 309 hs = (ctx->xmax-ctx->xmin)/2.0*(ctx->hratio+1.0)/Mx; 310 hf = (ctx->xmax-ctx->xmin)/2.0*(1.0+1.0/ctx->hratio)/Mx; 311 ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr); /* X is solution vector which does not contain ghost points */ 312 ierr = DMGlobalToLocalEnd (da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr); 313 ierr = VecZeroEntries(F);CHKERRQ(ierr); /* F is the right hand side function corresponds to center points */ 314 ierr = DMDAVecGetArray(da,Xloc,&x);CHKERRQ(ierr); 315 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 316 ierr = DMDAGetCorners(da,&xs,0,0,&xm,0,0);CHKERRQ(ierr); 317 ierr = PetscMalloc4(dof,&r,dof,&min,dof,&alpha,dof,&gamma);CHKERRQ(ierr); 318 319 if (ctx->bctype == FVBC_OUTFLOW) { 320 for (i=xs-2; i<0; i++) { 321 for (j=0; j<dof; j++) x[i*dof+j] = x[j]; 322 } 323 for (i=Mx; i<xs+xm+2; i++) { 324 for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j]; 325 } 326 } 327 328 for (i=xs; i<xs+xm+1; i++) { 329 PetscReal maxspeed; 330 PetscScalar *u; 331 if (i < sf) { 332 u = &ctx->u[0]; 333 alpha[0] = 1.0/6.0; 334 gamma[0] = 1.0/3.0; 335 for (j=0; j<dof; j++) { 336 r[j] = (x[i*dof+j] - x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]); 337 min[j] = PetscMin(r[j],2.0); 338 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]); 339 } 340 ierr = (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);CHKERRQ(ierr); 341 if (i > xs) { 342 for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hs; 343 } 344 if (i < xs+xm) { 345 for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hs; 346 islow++; 347 } 348 } else if (i == sf) { 349 u = &ctx->u[0]; 350 alpha[0] = hs*hf/(hs+hs)/(hs+hs+hf); 351 gamma[0] = hs*(hs+hs)/(hs+hf)/(hs+hs+hf); 352 for (j=0; j<dof; j++) { 353 r[j] = (x[i*dof+j] - x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]); 354 min[j] = PetscMin(r[j],2.0); 355 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]); 356 } 357 ierr = (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);CHKERRQ(ierr); 358 if (i > xs) { 359 for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hs; 360 } 361 } else if (i == fs) { 362 u = &ctx->u[0]; 363 alpha[0] = hf*hs/(hf+hf)/(hf+hf+hs); 364 gamma[0] = hf*(hf+hf)/(hf+hs)/(hf+hf+hs); 365 for (j=0; j<dof; j++) { 366 r[j] = (x[i*dof+j] - x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]); 367 min[j] = PetscMin(r[j],2.0); 368 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]); 369 } 370 ierr = (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);CHKERRQ(ierr); 371 if (i < xs+xm) { 372 for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hs; 373 islow++; 374 } 375 } else if (i == fs+1) { 376 u = &ctx->u[0]; 377 alpha[0] = hs*hs/(hf+hs)/(hf+hs+hs); 378 gamma[0] = hs*(hf+hs)/(hs+hs)/(hf+hs+hs); 379 for (j=0; j<dof; j++) { 380 r[j] = (x[i*dof+j] - x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]); 381 min[j] = PetscMin(r[j],2.0); 382 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]); 383 } 384 ierr = (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);CHKERRQ(ierr); 385 if (i > xs) { 386 for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hs; 387 } 388 if (i < xs+xm) { 389 for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hs; 390 islow++; 391 } 392 } else if (i > fs+1) { 393 u = &ctx->u[0]; 394 alpha[0] = 1.0/6.0; 395 gamma[0] = 1.0/3.0; 396 for (j=0; j<dof; j++) { 397 r[j] = (x[i*dof+j] - x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]); 398 min[j] = PetscMin(r[j],2.0); 399 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]); 400 } 401 ierr = (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);CHKERRQ(ierr); 402 if (i > xs) { 403 for (j=0; j<dof; j++) f[(islow-1)*dof+j] -= ctx->flux[j]/hs; 404 } 405 if (i < xs+xm) { 406 for (j=0; j<dof; j++) f[islow*dof+j] += ctx->flux[j]/hs; 407 islow++; 408 } 409 } 410 } 411 ierr = DMDAVecRestoreArray(da,Xloc,&x);CHKERRQ(ierr); 412 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 413 ierr = DMRestoreLocalVector(da,&Xloc);CHKERRQ(ierr); 414 ierr = PetscFree4(r,min,alpha,gamma);CHKERRQ(ierr); 415 PetscFunctionReturn(0); 416 } 417 418 static PetscErrorCode FVRHSFunctionfast(TS ts,PetscReal time,Vec X,Vec F,void *vctx) 419 { 420 FVCtx *ctx = (FVCtx*)vctx; 421 PetscErrorCode ierr; 422 PetscInt i,j,Mx,dof,xs,xm,ifast = 0,sf = ctx->sf,fs = ctx->fs; 423 PetscReal hf,hs; 424 PetscScalar *x,*f,*r,*min,*alpha,*gamma; 425 Vec Xloc; 426 DM da; 427 428 PetscFunctionBeginUser; 429 ierr = TSGetDM(ts,&da);CHKERRQ(ierr); 430 ierr = DMGetLocalVector(da,&Xloc);CHKERRQ(ierr); /* Xloc contains ghost points */ 431 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 */ 432 hs = (ctx->xmax-ctx->xmin)/2.0*(ctx->hratio+1.0)/Mx; 433 hf = (ctx->xmax-ctx->xmin)/2.0*(1.0+1.0/ctx->hratio)/Mx; 434 ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr); /* X is solution vector which does not contain ghost points */ 435 ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr); 436 ierr = VecZeroEntries(F);CHKERRQ(ierr); /* F is the right hand side function corresponds to center points */ 437 ierr = DMDAVecGetArray(da,Xloc,&x);CHKERRQ(ierr); 438 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 439 ierr = DMDAGetCorners(da,&xs,0,0,&xm,0,0);CHKERRQ(ierr); 440 ierr = PetscMalloc4(dof,&r,dof,&min,dof,&alpha,dof,&gamma);CHKERRQ(ierr); 441 442 if (ctx->bctype == FVBC_OUTFLOW) { 443 for (i=xs-2; i<0; i++) { 444 for (j=0; j<dof; j++) x[i*dof+j] = x[j]; 445 } 446 for (i=Mx; i<xs+xm+2; i++) { 447 for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j]; 448 } 449 } 450 451 for (i=xs; i<xs+xm+1; i++) { 452 PetscReal maxspeed; 453 PetscScalar *u; 454 if (i == sf) { 455 u = &ctx->u[0]; 456 alpha[0] = hs*hf/(hs+hs)/(hs+hs+hf); 457 gamma[0] = hs*(hs+hs)/(hs+hf)/(hs+hs+hf); 458 for (j=0; j<dof; j++) { 459 r[j] = (x[i*dof+j]-x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]); 460 min[j] = PetscMin(r[j],2.0); 461 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]); 462 } 463 ierr = (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);CHKERRQ(ierr); 464 if (i < xs+xm) { 465 for (j=0; j<dof; j++) f[ifast*dof+j] += ctx->flux[j]/hf; 466 ifast++; 467 } 468 } else if (i == sf+1) { 469 u = &ctx->u[0]; 470 alpha[0] = hf*hf/(hs+hf)/(hs+hf+hf); 471 gamma[0] = hf*(hs+hf)/(hf+hf)/(hs+hf+hf); 472 for (j=0; j<dof; j++) { 473 r[j] = (x[i*dof+j]-x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]); 474 min[j] = PetscMin(r[j],2.0); 475 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]); 476 } 477 ierr = (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);CHKERRQ(ierr); 478 if (i > xs) { 479 for (j=0; j<dof; j++) f[(ifast-1)*dof+j] -= ctx->flux[j]/hf; 480 } 481 if (i < xs+xm) { 482 for (j=0; j<dof; j++) f[ifast*dof+j] += ctx->flux[j]/hf; 483 ifast++; 484 } 485 } else if (i > sf+1 && i < fs) { 486 u = &ctx->u[0]; 487 alpha[0] = 1.0/6.0; 488 gamma[0] = 1.0/3.0; 489 for (j=0; j<dof; j++) { 490 r[j] = (x[i*dof+j]-x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]); 491 min[j] = PetscMin(r[j],2.0); 492 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]); 493 } 494 ierr = (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);CHKERRQ(ierr); 495 if (i > xs) { 496 for (j=0; j<dof; j++) f[(ifast-1)*dof+j] -= ctx->flux[j]/hf; 497 } 498 if (i < xs+xm) { 499 for (j=0; j<dof; j++) f[ifast*dof+j] += ctx->flux[j]/hf; 500 ifast++; 501 } 502 } else if (i == fs) { 503 u = &ctx->u[0]; 504 alpha[0] = hf*hs/(hf+hf)/(hf+hf+hs); 505 gamma[0] = hf*(hf+hf)/(hf+hs)/(hf+hf+hs); 506 for (j=0; j<dof; j++) { 507 r[j] = (x[i*dof+j]-x[(i-1)*dof+j])/(x[(i-1)*dof+j]-x[(i-2)*dof+j]); 508 min[j] = PetscMin(r[j],2.0); 509 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]); 510 } 511 ierr = (*ctx->physics.flux)(ctx->physics.user,u,ctx->flux,&maxspeed);CHKERRQ(ierr); 512 if (i > xs) { 513 for (j=0; j<dof; j++) f[(ifast-1)*dof+j] -= ctx->flux[j]/hf; 514 } 515 } 516 } 517 ierr = DMDAVecRestoreArray(da,Xloc,&x);CHKERRQ(ierr); 518 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 519 ierr = DMRestoreLocalVector(da,&Xloc);CHKERRQ(ierr); 520 ierr = PetscFree4(r,min,alpha,gamma);CHKERRQ(ierr); 521 PetscFunctionReturn(0); 522 } 523 524 /* --------------------------------- Finite Volume Solver for slow components ----------------------------------- */ 525 526 PetscErrorCode FVSample(FVCtx *ctx,DM da,PetscReal time,Vec U) 527 { 528 PetscErrorCode ierr; 529 PetscScalar *u,*uj,xj,xi; 530 PetscInt i,j,k,dof,xs,xm,Mx,count_slow,count_fast; 531 const PetscInt N=200; 532 533 PetscFunctionBeginUser; 534 if (!ctx->physics.sample) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Physics has not provided a sampling function"); 535 ierr = DMDAGetInfo(da,0,&Mx,0,0,0,0,0,&dof,0,0,0,0,0);CHKERRQ(ierr); 536 ierr = DMDAGetCorners(da,&xs,0,0,&xm,0,0);CHKERRQ(ierr); 537 ierr = DMDAVecGetArray(da,U,&u);CHKERRQ(ierr); 538 ierr = PetscMalloc1(dof,&uj);CHKERRQ(ierr); 539 const PetscReal hs = (ctx->xmax-ctx->xmin)/2.0*(ctx->hratio+1.0)/Mx; 540 const PetscReal hf = (ctx->xmax-ctx->xmin)/2.0*(1.0+1.0/ctx->hratio)/Mx; 541 count_slow = Mx/(1+ctx->hratio); 542 count_fast = Mx-count_slow; 543 for (i=xs; i<xs+xm; i++) { 544 if (i*hs+0.5*hs<(ctx->xmax-ctx->xmin)*0.25) { 545 xi = ctx->xmin+0.5*hs+i*hs; 546 /* Integrate over cell i using trapezoid rule with N points. */ 547 for (k=0; k<dof; k++) u[i*dof+k] = 0; 548 for (j=0; j<N+1; j++) { 549 xj = xi+hs*(j-N/2)/(PetscReal)N; 550 ierr = (*ctx->physics.sample)(ctx->physics.user,ctx->initial,ctx->bctype,ctx->xmin,ctx->xmax,time,xj,uj);CHKERRQ(ierr); 551 for (k=0; k<dof; k++) u[i*dof+k] += ((j==0 || j==N) ? 0.5 : 1.0)*uj[k]/N; 552 } 553 } else if ((ctx->xmax-ctx->xmin)*0.25+(i-count_slow/2)*hf+0.5*hf<(ctx->xmax-ctx->xmin)*0.75) { 554 xi = ctx->xmin+(ctx->xmax-ctx->xmin)*0.25+0.5*hf+(i-count_slow/2)*hf; 555 /* Integrate over cell i using trapezoid rule with N points. */ 556 for (k=0; k<dof; k++) u[i*dof+k] = 0; 557 for (j=0; j<N+1; j++) { 558 xj = xi+hf*(j-N/2)/(PetscReal)N; 559 ierr = (*ctx->physics.sample)(ctx->physics.user,ctx->initial,ctx->bctype,ctx->xmin,ctx->xmax,time,xj,uj);CHKERRQ(ierr); 560 for (k=0; k<dof; k++) u[i*dof+k] += ((j==0 || j==N) ? 0.5 : 1.0)*uj[k]/N; 561 } 562 } else { 563 xi = ctx->xmin+(ctx->xmax-ctx->xmin)*0.75+0.5*hs+(i-count_slow/2-count_fast)*hs; 564 /* Integrate over cell i using trapezoid rule with N points. */ 565 for (k=0; k<dof; k++) u[i*dof+k] = 0; 566 for (j=0; j<N+1; j++) { 567 xj = xi+hs*(j-N/2)/(PetscReal)N; 568 ierr = (*ctx->physics.sample)(ctx->physics.user,ctx->initial,ctx->bctype,ctx->xmin,ctx->xmax,time,xj,uj);CHKERRQ(ierr); 569 for (k=0; k<dof; k++) u[i*dof+k] += ((j==0 || j==N) ? 0.5 : 1.0)*uj[k]/N; 570 } 571 } 572 } 573 ierr = DMDAVecRestoreArray(da,U,&u);CHKERRQ(ierr); 574 ierr = PetscFree(uj);CHKERRQ(ierr); 575 PetscFunctionReturn(0); 576 } 577 578 static PetscErrorCode SolutionStatsView(DM da,Vec X,PetscViewer viewer) 579 { 580 PetscErrorCode ierr; 581 PetscReal xmin,xmax; 582 PetscScalar sum,tvsum,tvgsum; 583 const PetscScalar *x; 584 PetscInt imin,imax,Mx,i,j,xs,xm,dof; 585 Vec Xloc; 586 PetscBool iascii; 587 588 PetscFunctionBeginUser; 589 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 590 if (iascii) { 591 /* PETSc lacks a function to compute total variation norm (difficult in multiple dimensions), we do it here */ 592 ierr = DMGetLocalVector(da,&Xloc);CHKERRQ(ierr); 593 ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr); 594 ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr); 595 ierr = DMDAVecGetArrayRead(da,Xloc,(void*)&x);CHKERRQ(ierr); 596 ierr = DMDAGetCorners(da,&xs,0,0,&xm,0,0);CHKERRQ(ierr); 597 ierr = DMDAGetInfo(da,0,&Mx,0,0,0,0,0,&dof,0,0,0,0,0);CHKERRQ(ierr); 598 tvsum = 0; 599 for (i=xs; i<xs+xm; i++) { 600 for (j=0; j<dof; j++) tvsum += PetscAbsScalar(x[i*dof+j]-x[(i-1)*dof+j]); 601 } 602 ierr = MPI_Allreduce(&tvsum,&tvgsum,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)da));CHKERRMPI(ierr); 603 ierr = DMDAVecRestoreArrayRead(da,Xloc,(void*)&x);CHKERRQ(ierr); 604 ierr = DMRestoreLocalVector(da,&Xloc);CHKERRQ(ierr); 605 606 ierr = VecMin(X,&imin,&xmin);CHKERRQ(ierr); 607 ierr = VecMax(X,&imax,&xmax);CHKERRQ(ierr); 608 ierr = VecSum(X,&sum);CHKERRQ(ierr); 609 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); 610 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type not supported"); 611 PetscFunctionReturn(0); 612 } 613 614 static PetscErrorCode SolutionErrorNorms(FVCtx *ctx,DM da,PetscReal t,Vec X,PetscReal *nrm1) 615 { 616 PetscErrorCode ierr; 617 Vec Y; 618 PetscInt i,Mx,count_slow=0,count_fast=0; 619 const PetscScalar *ptr_X,*ptr_Y; 620 621 PetscFunctionBeginUser; 622 ierr = VecGetSize(X,&Mx);CHKERRQ(ierr); 623 ierr = VecDuplicate(X,&Y);CHKERRQ(ierr); 624 ierr = FVSample(ctx,da,t,Y);CHKERRQ(ierr); 625 const PetscReal hs = (ctx->xmax-ctx->xmin)/2.0*(ctx->hratio+1.0)/Mx; 626 const PetscReal hf = (ctx->xmax-ctx->xmin)/2.0*(1.0+1.0/ctx->hratio)/Mx; 627 count_slow = (PetscReal)Mx/(1.0+ctx->hratio); 628 count_fast = Mx-count_slow; 629 ierr = VecGetArrayRead(X,&ptr_X);CHKERRQ(ierr); 630 ierr = VecGetArrayRead(Y,&ptr_Y);CHKERRQ(ierr); 631 for (i=0; i<Mx; i++) { 632 if (i < count_slow/2 || i > count_slow/2+count_fast-1) *nrm1 += hs*PetscAbs(ptr_X[i]-ptr_Y[i]); 633 else *nrm1 += hf*PetscAbs(ptr_X[i]-ptr_Y[i]); 634 } 635 ierr = VecRestoreArrayRead(X,&ptr_X);CHKERRQ(ierr); 636 ierr = VecRestoreArrayRead(Y,&ptr_Y);CHKERRQ(ierr); 637 ierr = VecDestroy(&Y);CHKERRQ(ierr); 638 PetscFunctionReturn(0); 639 } 640 641 int main(int argc,char *argv[]) 642 { 643 char physname[256] = "advect",final_fname[256] = "solution.m"; 644 PetscFunctionList physics = 0; 645 MPI_Comm comm; 646 TS ts; 647 DM da; 648 Vec X,X0,R; 649 FVCtx ctx; 650 PetscInt i,k,dof,xs,xm,Mx,draw = 0,count_slow,count_fast,islow = 0,ifast = 0,*index_slow,*index_fast; 651 PetscBool view_final = PETSC_FALSE; 652 PetscReal ptime; 653 PetscErrorCode ierr; 654 655 ierr = PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr; 656 comm = PETSC_COMM_WORLD; 657 ierr = PetscMemzero(&ctx,sizeof(ctx));CHKERRQ(ierr); 658 659 /* Register physical models to be available on the command line */ 660 ierr = PetscFunctionListAdd(&physics,"advect",PhysicsCreate_Advect);CHKERRQ(ierr); 661 662 ctx.comm = comm; 663 ctx.cfl = 0.9; 664 ctx.bctype = FVBC_PERIODIC; 665 ctx.xmin = -1.0; 666 ctx.xmax = 1.0; 667 ierr = PetscOptionsBegin(comm,NULL,"Finite Volume solver options","");CHKERRQ(ierr); 668 ierr = PetscOptionsReal("-xmin","X min","",ctx.xmin,&ctx.xmin,NULL);CHKERRQ(ierr); 669 ierr = PetscOptionsReal("-xmax","X max","",ctx.xmax,&ctx.xmax,NULL);CHKERRQ(ierr); 670 ierr = PetscOptionsInt("-draw","Draw solution vector, bitwise OR of (1=initial,2=final,4=final error)","",draw,&draw,NULL);CHKERRQ(ierr); 671 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); 672 ierr = PetscOptionsInt("-initial","Initial condition (depends on the physics)","",ctx.initial,&ctx.initial,NULL);CHKERRQ(ierr); 673 ierr = PetscOptionsBool("-exact","Compare errors with exact solution","",ctx.exact,&ctx.exact,NULL);CHKERRQ(ierr); 674 ierr = PetscOptionsBool("-simulation","Compare errors with reference solution","",ctx.simulation,&ctx.simulation,NULL);CHKERRQ(ierr); 675 ierr = PetscOptionsReal("-cfl","CFL number to time step at","",ctx.cfl,&ctx.cfl,NULL);CHKERRQ(ierr); 676 ierr = PetscOptionsEnum("-bc_type","Boundary condition","",FVBCTypes,(PetscEnum)ctx.bctype,(PetscEnum*)&ctx.bctype,NULL);CHKERRQ(ierr); 677 ierr = PetscOptionsInt("-hratio","Spacing ratio","",ctx.hratio,&ctx.hratio,NULL);CHKERRQ(ierr); 678 ierr = PetscOptionsEnd();CHKERRQ(ierr); 679 680 /* Choose the physics from the list of registered models */ 681 { 682 PetscErrorCode (*r)(FVCtx*); 683 ierr = PetscFunctionListFind(physics,physname,&r);CHKERRQ(ierr); 684 if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Physics '%s' not found",physname); 685 /* Create the physics, will set the number of fields and their names */ 686 ierr = (*r)(&ctx);CHKERRQ(ierr); 687 } 688 689 /* Create a DMDA to manage the parallel grid */ 690 ierr = DMDACreate1d(comm,DM_BOUNDARY_PERIODIC,50,ctx.physics.dof,2,NULL,&da);CHKERRQ(ierr); 691 ierr = DMSetFromOptions(da);CHKERRQ(ierr); 692 ierr = DMSetUp(da);CHKERRQ(ierr); 693 /* Inform the DMDA of the field names provided by the physics. */ 694 /* The names will be shown in the title bars when run with -ts_monitor_draw_solution */ 695 for (i=0; i<ctx.physics.dof; i++) { 696 ierr = DMDASetFieldName(da,i,ctx.physics.fieldname[i]);CHKERRQ(ierr); 697 } 698 ierr = DMDAGetInfo(da,0,&Mx,0,0,0,0,0,&dof,0,0,0,0,0);CHKERRQ(ierr); 699 ierr = DMDAGetCorners(da,&xs,0,0,&xm,0,0);CHKERRQ(ierr); 700 701 /* Set coordinates of cell centers */ 702 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); 703 704 /* Allocate work space for the Finite Volume solver (so it doesn't have to be reallocated on each function evaluation) */ 705 ierr = PetscMalloc3(dof,&ctx.u,dof,&ctx.flux,dof,&ctx.speeds);CHKERRQ(ierr); 706 707 /* Create a vector to store the solution and to save the initial state */ 708 ierr = DMCreateGlobalVector(da,&X);CHKERRQ(ierr); 709 ierr = VecDuplicate(X,&X0);CHKERRQ(ierr); 710 ierr = VecDuplicate(X,&R);CHKERRQ(ierr); 711 712 /* create index for slow parts and fast parts*/ 713 count_slow = Mx/(1+ctx.hratio); 714 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"); 715 count_fast = Mx-count_slow; 716 ctx.sf = count_slow/2; 717 ctx.fs = ctx.sf + count_fast; 718 ierr = PetscMalloc1(xm*dof,&index_slow);CHKERRQ(ierr); 719 ierr = PetscMalloc1(xm*dof,&index_fast);CHKERRQ(ierr); 720 for (i=xs; i<xs+xm; i++) { 721 if (i < count_slow/2 || i > count_slow/2+count_fast-1) 722 for (k=0; k<dof; k++) index_slow[islow++] = i*dof+k; 723 else 724 for (k=0; k<dof; k++) index_fast[ifast++] = i*dof+k; 725 } 726 ierr = ISCreateGeneral(PETSC_COMM_WORLD,islow,index_slow,PETSC_COPY_VALUES,&ctx.iss);CHKERRQ(ierr); 727 ierr = ISCreateGeneral(PETSC_COMM_WORLD,ifast,index_fast,PETSC_COPY_VALUES,&ctx.isf);CHKERRQ(ierr); 728 729 /* Create a time-stepping object */ 730 ierr = TSCreate(comm,&ts);CHKERRQ(ierr); 731 ierr = TSSetDM(ts,da);CHKERRQ(ierr); 732 ierr = TSSetRHSFunction(ts,R,FVRHSFunction,&ctx);CHKERRQ(ierr); 733 ierr = TSRHSSplitSetIS(ts,"slow",ctx.iss);CHKERRQ(ierr); 734 ierr = TSRHSSplitSetIS(ts,"fast",ctx.isf);CHKERRQ(ierr); 735 ierr = TSRHSSplitSetRHSFunction(ts,"slow",NULL,FVRHSFunctionslow,&ctx);CHKERRQ(ierr); 736 ierr = TSRHSSplitSetRHSFunction(ts,"fast",NULL,FVRHSFunctionfast,&ctx);CHKERRQ(ierr); 737 738 ierr = TSSetType(ts,TSMPRK);CHKERRQ(ierr); 739 ierr = TSSetMaxTime(ts,10);CHKERRQ(ierr); 740 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr); 741 742 /* Compute initial conditions and starting time step */ 743 ierr = FVSample(&ctx,da,0,X0);CHKERRQ(ierr); 744 ierr = FVRHSFunction(ts,0,X0,X,(void*)&ctx);CHKERRQ(ierr); /* Initial function evaluation, only used to determine max speed */ 745 ierr = VecCopy(X0,X);CHKERRQ(ierr); /* The function value was not used so we set X=X0 again */ 746 ierr = TSSetTimeStep(ts,ctx.cfl/ctx.cfl_idt);CHKERRQ(ierr); 747 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); /* Take runtime options */ 748 ierr = SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 749 { 750 PetscInt steps; 751 PetscScalar mass_initial,mass_final,mass_difference,mass_differenceg; 752 const PetscScalar *ptr_X,*ptr_X0; 753 const PetscReal hs = (ctx.xmax-ctx.xmin)/2.0/count_slow; 754 const PetscReal hf = (ctx.xmax-ctx.xmin)/2.0/count_fast; 755 ierr = TSSolve(ts,X);CHKERRQ(ierr); 756 ierr = TSGetSolveTime(ts,&ptime);CHKERRQ(ierr); 757 ierr = TSGetStepNumber(ts,&steps);CHKERRQ(ierr); 758 /* calculate the total mass at initial time and final time */ 759 mass_initial = 0.0; 760 mass_final = 0.0; 761 ierr = DMDAVecGetArrayRead(da,X0,(void*)&ptr_X0);CHKERRQ(ierr); 762 ierr = DMDAVecGetArrayRead(da,X,(void*)&ptr_X);CHKERRQ(ierr); 763 for (i=xs; i<xs+xm; i++) { 764 if (i < ctx.sf || i > ctx.fs-1) { 765 for (k=0; k<dof; k++) { 766 mass_initial = mass_initial+hs*ptr_X0[i*dof+k]; 767 mass_final = mass_final+hs*ptr_X[i*dof+k]; 768 } 769 } else { 770 for (k=0; k<dof; k++) { 771 mass_initial = mass_initial+hf*ptr_X0[i*dof+k]; 772 mass_final = mass_final+hf*ptr_X[i*dof+k]; 773 } 774 } 775 } 776 ierr = DMDAVecRestoreArrayRead(da,X0,(void*)&ptr_X0);CHKERRQ(ierr); 777 ierr = DMDAVecRestoreArrayRead(da,X,(void*)&ptr_X);CHKERRQ(ierr); 778 mass_difference = mass_final-mass_initial; 779 ierr = MPI_Allreduce(&mass_difference,&mass_differenceg,1,MPIU_SCALAR,MPIU_SUM,comm);CHKERRMPI(ierr); 780 ierr = PetscPrintf(comm,"Mass difference %g\n",(double)mass_differenceg);CHKERRQ(ierr); 781 ierr = PetscPrintf(comm,"Final time %g, steps %D\n",(double)ptime,steps);CHKERRQ(ierr); 782 if (ctx.exact) { 783 PetscReal nrm1 = 0; 784 ierr = SolutionErrorNorms(&ctx,da,ptime,X,&nrm1);CHKERRQ(ierr); 785 ierr = PetscPrintf(comm,"Error ||x-x_e||_1 %g\n",(double)nrm1);CHKERRQ(ierr); 786 } 787 if (ctx.simulation) { 788 PetscReal nrm1 = 0; 789 PetscViewer fd; 790 char filename[PETSC_MAX_PATH_LEN] = "binaryoutput"; 791 Vec XR; 792 PetscBool flg; 793 const PetscScalar *ptr_XR; 794 ierr = PetscOptionsGetString(NULL,NULL,"-f",filename,sizeof(filename),&flg);CHKERRQ(ierr); 795 if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Must indicate binary file with the -f option"); 796 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&fd);CHKERRQ(ierr); 797 ierr = VecDuplicate(X0,&XR);CHKERRQ(ierr); 798 ierr = VecLoad(XR,fd);CHKERRQ(ierr); 799 ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr); 800 ierr = VecGetArrayRead(X,&ptr_X);CHKERRQ(ierr); 801 ierr = VecGetArrayRead(XR,&ptr_XR);CHKERRQ(ierr); 802 for (i=0; i<Mx; i++) { 803 if (i < count_slow/2 || i > count_slow/2+count_fast-1) nrm1 = nrm1 + hs*PetscAbs(ptr_X[i]-ptr_XR[i]); 804 else nrm1 = nrm1 + hf*PetscAbs(ptr_X[i]-ptr_XR[i]); 805 } 806 ierr = VecRestoreArrayRead(X,&ptr_X);CHKERRQ(ierr); 807 ierr = VecRestoreArrayRead(XR,&ptr_XR);CHKERRQ(ierr); 808 ierr = PetscPrintf(comm,"Error ||x-x_e||_1 %g\n",(double)nrm1);CHKERRQ(ierr); 809 ierr = VecDestroy(&XR);CHKERRQ(ierr); 810 } 811 } 812 813 ierr = SolutionStatsView(da,X,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 814 if (draw & 0x1) { ierr = VecView(X0,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr); } 815 if (draw & 0x2) { ierr = VecView(X,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr); } 816 if (draw & 0x4) { 817 Vec Y; 818 ierr = VecDuplicate(X,&Y);CHKERRQ(ierr); 819 ierr = FVSample(&ctx,da,ptime,Y);CHKERRQ(ierr); 820 ierr = VecAYPX(Y,-1,X);CHKERRQ(ierr); 821 ierr = VecView(Y,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr); 822 ierr = VecDestroy(&Y);CHKERRQ(ierr); 823 } 824 825 if (view_final) { 826 PetscViewer viewer; 827 ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD,final_fname,&viewer);CHKERRQ(ierr); 828 ierr = PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);CHKERRQ(ierr); 829 ierr = VecView(X,viewer);CHKERRQ(ierr); 830 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 831 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 832 } 833 834 /* Clean up */ 835 ierr = (*ctx.physics.destroy)(ctx.physics.user);CHKERRQ(ierr); 836 for (i=0; i<ctx.physics.dof; i++) {ierr = PetscFree(ctx.physics.fieldname[i]);CHKERRQ(ierr);} 837 ierr = PetscFree3(ctx.u,ctx.flux,ctx.speeds);CHKERRQ(ierr); 838 ierr = ISDestroy(&ctx.iss);CHKERRQ(ierr); 839 ierr = ISDestroy(&ctx.isf);CHKERRQ(ierr); 840 ierr = VecDestroy(&X);CHKERRQ(ierr); 841 ierr = VecDestroy(&X0);CHKERRQ(ierr); 842 ierr = VecDestroy(&R);CHKERRQ(ierr); 843 ierr = DMDestroy(&da);CHKERRQ(ierr); 844 ierr = TSDestroy(&ts);CHKERRQ(ierr); 845 ierr = PetscFree(index_slow);CHKERRQ(ierr); 846 ierr = PetscFree(index_fast);CHKERRQ(ierr); 847 ierr = PetscFunctionListDestroy(&physics);CHKERRQ(ierr); 848 ierr = PetscFinalize(); 849 return ierr; 850 } 851 852 /*TEST 853 854 build: 855 requires: !complex 856 857 test: 858 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 859 860 test: 861 suffix: 2 862 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 863 output_file: output/ex7_1.out 864 865 test: 866 suffix: 3 867 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 868 869 test: 870 suffix: 4 871 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 872 output_file: output/ex7_3.out 873 874 test: 875 suffix: 5 876 nsize: 2 877 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 878 output_file: output/ex7_3.out 879 TEST*/ 880