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