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