1 #include "finitevolume1d.h" 2 #include <petscdmda.h> 3 #include <petscdraw.h> 4 #include <petsc/private/tsimpl.h> 5 6 #include <petsc/private/kernels/blockinvert.h> /* For the Kernel_*_gets_* stuff for BAIJ */ 7 const char *FVBCTypes[] = {"PERIODIC","OUTFLOW","FVBCType","FVBC_",0}; 8 9 PETSC_STATIC_INLINE PetscReal Sgn(PetscReal a) { return (a<0) ? -1 : 1; } 10 PETSC_STATIC_INLINE PetscReal Abs(PetscReal a) { return (a<0) ? 0 : a; } 11 PETSC_STATIC_INLINE PetscReal Sqr(PetscReal a) { return a*a; } 12 //PETSC_STATIC_INLINE PetscReal MaxAbs(PetscReal a,PetscReal b) { return (PetscAbs(a) > PetscAbs(b)) ? a : b; } 13 PETSC_UNUSED PETSC_STATIC_INLINE PetscReal MinAbs(PetscReal a,PetscReal b) { return (PetscAbs(a) < PetscAbs(b)) ? a : b; } 14 PETSC_STATIC_INLINE PetscReal MinMod2(PetscReal a,PetscReal b) { return (a*b<0) ? 0 : Sgn(a)*PetscMin(PetscAbs(a),PetscAbs(b)); } 15 PETSC_STATIC_INLINE PetscReal MaxMod2(PetscReal a,PetscReal b) { return (a*b<0) ? 0 : Sgn(a)*PetscMax(PetscAbs(a),PetscAbs(b)); } 16 PETSC_STATIC_INLINE PetscReal MinMod3(PetscReal a,PetscReal b,PetscReal c) {return (a*b<0 || a*c<0) ? 0 : Sgn(a)*PetscMin(PetscAbs(a),PetscMin(PetscAbs(b),PetscAbs(c))); } 17 18 //PETSC_STATIC_INLINE PetscReal RangeMod(PetscReal a,PetscReal xmin,PetscReal xmax) { PetscReal range = xmax-xmin; return xmin +PetscFmodReal(range+PetscFmodReal(a,range),range); } 19 20 21 /* ----------------------- Lots of limiters, these could go in a separate library ------------------------- */ 22 void Limit_Upwind(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) 23 { 24 PetscInt i; 25 for (i=0; i<info->m; i++) lmt[i] = 0; 26 } 27 void Limit_LaxWendroff(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) 28 { 29 PetscInt i; 30 for (i=0; i<info->m; i++) lmt[i] = jR[i]; 31 } 32 void Limit_BeamWarming(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) 33 { 34 PetscInt i; 35 for (i=0; i<info->m; i++) lmt[i] = jL[i]; 36 } 37 void Limit_Fromm(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) 38 { 39 PetscInt i; 40 for (i=0; i<info->m; i++) lmt[i] = 0.5*(jL[i]+jR[i]); 41 } 42 void Limit_Minmod(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) 43 { 44 PetscInt i; 45 for (i=0; i<info->m; i++) lmt[i] = MinMod2(jL[i],jR[i]); 46 } 47 void Limit_Superbee(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) 48 { 49 PetscInt i; 50 for (i=0; i<info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i],2*jR[i]),MinMod2(2*jL[i],jR[i])); 51 } 52 void Limit_MC(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) 53 { 54 PetscInt i; 55 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i],0.5*(jL[i]+jR[i]),2*jR[i]); 56 } 57 void Limit_VanLeer(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) 58 { /* phi = (t + abs(t)) / (1 + abs(t)) */ 59 PetscInt i; 60 for (i=0; i<info->m; i++) lmt[i] = (jL[i]*Abs(jR[i])+Abs(jL[i])*jR[i])/(Abs(jL[i])+Abs(jR[i])+1e-15); 61 } 62 void Limit_VanAlbada(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) /* differentiable */ 63 { /* phi = (t + t^2) / (1 + t^2) */ 64 PetscInt i; 65 for (i=0; i<info->m; i++) lmt[i] = (jL[i]*Sqr(jR[i])+Sqr(jL[i])*jR[i])/(Sqr(jL[i])+Sqr(jR[i])+1e-15); 66 } 67 void Limit_VanAlbadaTVD(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) 68 { /* phi = (t + t^2) / (1 + t^2) */ 69 PetscInt i; 70 for (i=0; i<info->m; i++) lmt[i] = (jL[i]*jR[i]<0) ? 0 : (jL[i]*Sqr(jR[i])+Sqr(jL[i])*jR[i])/(Sqr(jL[i])+Sqr(jR[i])+1e-15); 71 } 72 void Limit_Koren(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) /* differentiable */ 73 { /* phi = (t + 2*t^2) / (2 - t + 2*t^2) */ 74 PetscInt i; 75 for (i=0; i<info->m; i++) lmt[i] = ((jL[i]*Sqr(jR[i])+2*Sqr(jL[i])*jR[i])/(2*Sqr(jL[i])-jL[i]*jR[i]+2*Sqr(jR[i])+1e-15)); 76 } 77 void Limit_KorenSym(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) /* differentiable */ 78 { /* Symmetric version of above */ 79 PetscInt i; 80 for (i=0; i<info->m; i++) lmt[i] = (1.5*(jL[i]*Sqr(jR[i])+Sqr(jL[i])*jR[i])/(2*Sqr(jL[i])-jL[i]*jR[i]+2*Sqr(jR[i])+1e-15)); 81 } 82 void Limit_Koren3(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) 83 { /* Eq 11 of Cada-Torrilhon 2009 */ 84 PetscInt i; 85 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i],(jL[i]+2*jR[i])/3,2*jR[i]); 86 } 87 static PetscReal CadaTorrilhonPhiHatR_Eq13(PetscReal L,PetscReal R) 88 { 89 return PetscMax(0,PetscMin((L+2*R)/3,PetscMax(-0.5*L,PetscMin(2*L,PetscMin((L+2*R)/3,1.6*R))))); 90 } 91 void Limit_CadaTorrilhon2(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) 92 { /* Cada-Torrilhon 2009, Eq 13 */ 93 PetscInt i; 94 for (i=0; i<info->m; i++) lmt[i] = CadaTorrilhonPhiHatR_Eq13(jL[i],jR[i]); 95 } 96 void Limit_CadaTorrilhon3R(PetscReal r,LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) 97 { /* Cada-Torrilhon 2009, Eq 22 */ 98 /* They recommend 0.001 < r < 1, but larger values are more accurate in smooth regions */ 99 const PetscReal eps = 1e-7,hx = info->hx; 100 PetscInt i; 101 for (i=0; i<info->m; i++) { 102 const PetscReal eta = (Sqr(jL[i])+Sqr(jR[i]))/Sqr(r*hx); 103 lmt[i] = ((eta < 1-eps) ? (jL[i]+2*jR[i])/3 : ((eta>1+eps) ? CadaTorrilhonPhiHatR_Eq13(jL[i],jR[i]) : 0.5*((1-(eta-1)/eps)*(jL[i]+2*jR[i])/3+(1+(eta+1)/eps)*CadaTorrilhonPhiHatR_Eq13(jL[i],jR[i])))); 104 } 105 } 106 void Limit_CadaTorrilhon3R0p1(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) 107 { 108 Limit_CadaTorrilhon3R(0.1,info,jL,jR,lmt); 109 } 110 void Limit_CadaTorrilhon3R1(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) 111 { 112 Limit_CadaTorrilhon3R(1,info,jL,jR,lmt); 113 } 114 void Limit_CadaTorrilhon3R10(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) 115 { 116 Limit_CadaTorrilhon3R(10,info,jL,jR,lmt); 117 } 118 void Limit_CadaTorrilhon3R100(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) 119 { 120 Limit_CadaTorrilhon3R(100,info,jL,jR,lmt); 121 } 122 123 /* ----------------------- Limiters for split systems ------------------------- */ 124 125 void Limit2_Upwind(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,const PetscInt sf,const PetscInt fs,PetscInt n,PetscScalar *lmt) 126 { 127 PetscInt i; 128 for (i=0; i<info->m; i++) lmt[i] = 0; 129 } 130 void Limit2_LaxWendroff(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,const PetscInt sf,const PetscInt fs,PetscInt n,PetscScalar *lmt) 131 { 132 PetscInt i; 133 if (n < sf-1) { /* slow components */ 134 for (i=0; i<info->m; i++) lmt[i] = jR[i]/info->hxs; 135 } else if (n == sf-1) { /* slow component which is next to fast components */ 136 for (i=0; i<info->m; i++) lmt[i] = jR[i]/(info->hxs/2.0+info->hxf/2.0); 137 } else if (n == sf) { /* fast component which is next to slow components */ 138 for (i=0; i<info->m; i++) lmt[i] = jR[i]/info->hxf; 139 } else if (n > sf && n < fs-1) { /* fast components */ 140 for (i=0; i<info->m; i++) lmt[i] = jR[i]/info->hxf; 141 } else if (n == fs-1) { /* fast component next to slow components */ 142 for (i=0; i<info->m; i++) lmt[i] = jR[i]/(info->hxf/2.0+info->hxs/2.0); 143 } else if (n == fs) { /* slow component next to fast components */ 144 for (i=0; i<info->m; i++) lmt[i] = jR[i]/info->hxs; 145 } else { /* slow components */ 146 for (i=0; i<info->m; i++) lmt[i] = jR[i]/info->hxs; 147 } 148 } 149 void Limit2_BeamWarming(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,const PetscInt sf,const PetscInt fs,PetscInt n,PetscScalar *lmt) 150 { 151 PetscInt i; 152 if (n < sf-1) { 153 for (i=0; i<info->m; i++) lmt[i] = jL[i]/info->hxs; 154 } else if (n == sf-1) { 155 for (i=0; i<info->m; i++) lmt[i] = jL[i]/info->hxs; 156 } else if (n == sf) { 157 for (i=0; i<info->m; i++) lmt[i] = jL[i]/(info->hxs/2.0+info->hxf/2.0); 158 } else if (n > sf && n < fs-1) { 159 for (i=0; i<info->m; i++) lmt[i] = jL[i]/info->hxf; 160 } else if (n == fs-1) { 161 for (i=0; i<info->m; i++) lmt[i] = jL[i]/info->hxf; 162 } else if (n == fs) { 163 for (i=0; i<info->m; i++) lmt[i] = jL[i]/(info->hxf/2.0+info->hxs/2.0); 164 } else { 165 for (i=0; i<info->m; i++) lmt[i] = jL[i]/info->hxs; 166 } 167 } 168 void Limit2_Fromm(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,const PetscInt sf,const PetscInt fs,PetscInt n,PetscScalar *lmt) 169 { 170 PetscInt i; 171 if (n < sf-1) { 172 for (i=0; i<info->m; i++) lmt[i] = 0.5*(jL[i]+jR[i])/info->hxs; 173 } else if (n == sf-1) { 174 for (i=0; i<info->m; i++) lmt[i] = 0.5*(jL[i]/info->hxs+jR[i]/(info->hxs/2.0+info->hxf/2.0)); 175 } else if (n == sf) { 176 for (i=0; i<info->m; i++) lmt[i] = 0.5*(jL[i]/(info->hxs/2.0+info->hxf/2.0)+jR[i]/info->hxf); 177 } else if (n > sf && n < fs-1) { 178 for (i=0; i<info->m; i++) lmt[i] = 0.5*(jL[i]+jR[i])/info->hxf; 179 } else if (n == fs-1) { 180 for (i=0; i<info->m; i++) lmt[i] = 0.5*(jL[i]/info->hxf+jR[i]/(info->hxf/2.0+info->hxs/2.0)); 181 } else if (n == fs) { 182 for (i=0; i<info->m; i++) lmt[i] = 0.5*(jL[i]/(info->hxf/2.0+info->hxs/2.0)+jR[i]/info->hxs); 183 } else { 184 for (i=0; i<info->m; i++) lmt[i] = 0.5*(jL[i]+jR[i])/info->hxs; 185 } 186 } 187 void Limit2_Minmod(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,const PetscInt sf,const PetscInt fs,PetscInt n,PetscScalar *lmt) 188 { 189 PetscInt i; 190 if (n < sf-1) { 191 for (i=0; i<info->m; i++) lmt[i] = MinMod2(jL[i],jR[i])/info->hxs; 192 } else if (n == sf-1) { 193 for (i=0; i<info->m; i++) lmt[i] = MinMod2(jL[i]/info->hxs,jR[i]/(info->hxs/2.0+info->hxf/2.0)); 194 } else if (n == sf) { 195 for (i=0; i<info->m; i++) lmt[i] = MinMod2(jL[i]/(info->hxs/2.0+info->hxf/2.0),jR[i]/info->hxf); 196 } else if (n > sf && n < fs-1) { 197 for (i=0; i<info->m; i++) lmt[i] = MinMod2(jL[i],jR[i])/info->hxf; 198 } else if (n == fs-1) { 199 for (i=0; i<info->m; i++) lmt[i] = MinMod2(jL[i]/info->hxf,jR[i]/(info->hxf/2.0+info->hxs/2.0)); 200 } else if (n == fs) { 201 for (i=0; i<info->m; i++) lmt[i] = MinMod2(jL[i]/(info->hxf/2.0+info->hxs/2.0),jR[i]/info->hxs); 202 } else { 203 for (i=0; i<info->m; i++) lmt[i] = MinMod2(jL[i],jR[i])/info->hxs; 204 } 205 } 206 void Limit2_Superbee(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,const PetscInt sf,const PetscInt fs,PetscInt n,PetscScalar *lmt) 207 { 208 PetscInt i; 209 if (n < sf-1) { 210 for (i=0; i<info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i],2*jR[i]),MinMod2(2*jL[i],jR[i]))/info->hxs; 211 } else if (n == sf-1) { 212 for (i=0; i<info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i]/info->hxs,2*jR[i]/(info->hxs/2.0+info->hxf/2.0)),MinMod2(2*jL[i]/info->hxs,jR[i]/(info->hxs/2.0+info->hxf/2.0))); 213 } else if (n == sf) { 214 for (i=0; i<info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i]/(info->hxs/2.0+info->hxf/2.0),2*jR[i]/info->hxf),MinMod2(2*jL[i]/(info->hxs/2.0+info->hxf/2.0),jR[i]/info->hxf)); 215 } else if (n > sf && n < fs-1) { 216 for (i=0; i<info->m; i++) lmt[i] = MinMod2(MinMod2(jL[i],2*jR[i]),MinMod2(2*jL[i],jR[i]))/info->hxf; 217 } else if (n == fs-1) { 218 for (i=0; i<info->m; i++) lmt[i] = MinMod2(MinMod2(jL[i]/info->hxf,2*jR[i]/(info->hxf/2.0+info->hxs/2.0)),MinMod2(2*jL[i]/info->hxf,jR[i]/(info->hxf/2.0+info->hxs/2.0))); 219 } else if (n == fs) { 220 for (i=0; i<info->m; i++) lmt[i] = MinMod2(MinMod2(jL[i]/(info->hxf/2.0+info->hxs/2.0),2*jR[i]/info->hxs),MinMod2(2*jL[i]/(info->hxf/2.0+info->hxs/2.0),jR[i]/info->hxs)); 221 } else { 222 for (i=0; i<info->m; i++) lmt[i] = MinMod2(MinMod2(jL[i],2*jR[i]),MinMod2(2*jL[i],jR[i]))/info->hxs; 223 } 224 } 225 void Limit2_MC(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,const PetscInt sf,const PetscInt fs,PetscInt n,PetscScalar *lmt) 226 { 227 PetscInt i; 228 if (n < sf-1) { 229 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i],0.5*(jL[i]+jR[i]),2*jR[i])/info->hxs; 230 } else if (n == sf-1) { 231 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i]/info->hxs,0.5*(jL[i]/info->hxs+jR[i]/(info->hxf/2.0+info->hxs/2.0)),2*jR[i]/(info->hxf/2.0+info->hxs/2.0)); 232 } else if (n == sf) { 233 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i]/(info->hxs/2.0+info->hxf/2.0),0.5*(jL[i]/(info->hxs/2.0+info->hxf/2.0)+jR[i]/info->hxf),2*jR[i]/info->hxf); 234 } else if (n > sf && n < fs-1) { 235 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i],0.5*(jL[i]+jR[i]),2*jR[i])/info->hxf; 236 } else if (n == fs-1) { 237 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i]/info->hxf,0.5*(jL[i]/info->hxf+jR[i]/(info->hxf/2.0+info->hxs/2.0)),2*jR[i]/(info->hxf/2.0+info->hxs/2.0)); 238 } else if (n == fs) { 239 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i]/(info->hxf/2.0+info->hxs/2.0),0.5*(jL[i]/(info->hxf/2.0+info->hxs/2.0)+jR[i]/info->hxs),2*jR[i]/info->hxs); 240 } else { 241 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i],0.5*(jL[i]+jR[i]),2*jR[i])/info->hxs; 242 } 243 } 244 void Limit2_Koren3(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,const PetscInt sf,const PetscInt fs,PetscInt n,PetscScalar *lmt) 245 { /* Eq 11 of Cada-Torrilhon 2009 */ 246 PetscInt i; 247 if (n < sf-1) { 248 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i],(jL[i]+2*jR[i])/3,2*jR[i])/info->hxs; 249 } else if (n == sf-1) { 250 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i]/info->hxs,(jL[i]/info->hxs+2*jR[i]/(info->hxf/2.0+info->hxs/2.0))/3,2*jR[i]/(info->hxf/2.0+info->hxs/2.0)); 251 } else if (n == sf) { 252 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i]/(info->hxs/2.0+info->hxf/2.0),(jL[i]/(info->hxs/2.0+info->hxf/2.0)+2*jR[i]/info->hxf)/3,2*jR[i]/info->hxf); 253 } else if (n > sf && n < fs-1) { 254 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i],(jL[i]+2*jR[i])/3,2*jR[i])/info->hxf; 255 } else if (n == fs-1) { 256 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i]/info->hxf,(jL[i]/info->hxf+2*jR[i]/(info->hxf/2.0+info->hxs/2.0))/3,2*jR[i]/(info->hxf/2.0+info->hxs/2.0)); 257 } else if (n == fs) { 258 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i]/(info->hxf/2.0+info->hxs/2.0),(jL[i]/(info->hxf/2.0+info->hxs/2.0)+2*jR[i]/info->hxs)/3,2*jR[i]/info->hxs); 259 } else { 260 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i],(jL[i]+2*jR[i])/3,2*jR[i])/info->hxs; 261 } 262 } 263 264 /* ---- Three-way splitting ---- */ 265 void Limit3_Upwind(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,const PetscInt sm,const PetscInt mf,const PetscInt fm,const PetscInt ms,PetscInt n,PetscScalar *lmt) 266 { 267 PetscInt i; 268 for (i=0; i<info->m; i++) lmt[i] = 0; 269 } 270 void Limit3_LaxWendroff(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,const PetscInt sm,const PetscInt mf,const PetscInt fm,const PetscInt ms,PetscInt n,PetscScalar *lmt) 271 { 272 PetscInt i; 273 if (n < sm-1 || n > ms) { /* slow components */ 274 for (i=0; i<info->m; i++) lmt[i] = jR[i]/info->hxs; 275 } else if (n == sm-1 || n == ms-1) { /* slow component which is next to medium components */ 276 for (i=0; i<info->m; i++) lmt[i] = jR[i]/(info->hxs/2.0+info->hxm/2.0); 277 } else if (n < mf-1 || n > fm) { /* medium components */ 278 for (i=0; i<info->m; i++) lmt[i] = jR[i]/info->hxm; 279 } else if (n == mf-1 || n == fm) { /* medium component next to fast components */ 280 for (i=0; i<info->m; i++) lmt[i] = jR[i]/(info->hxm/2.0+info->hxf/2.0); 281 } else { /* fast components */ 282 for (i=0; i<info->m; i++) lmt[i] = jR[i]/info->hxf; 283 } 284 } 285 void Limit3_BeamWarming(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,const PetscInt sm,const PetscInt mf,const PetscInt fm,const PetscInt ms,PetscInt n,PetscScalar *lmt) 286 { 287 PetscInt i; 288 if (n < sm || n > ms) { 289 for (i=0; i<info->m; i++) lmt[i] = jL[i]/info->hxs; 290 } else if (n == sm || n == ms) { 291 for (i=0; i<info->m; i++) lmt[i] = jL[i]/(info->hxs/2.0+info->hxf/2.0); 292 }else if (n < mf || n > fm) { 293 for (i=0; i<info->m; i++) lmt[i] = jL[i]/info->hxm; 294 } else if (n == mf || n == fm) { 295 for (i=0; i<info->m; i++) lmt[i] = jL[i]/(info->hxm/2.0+info->hxf/2.0); 296 } else { 297 for (i=0; i<info->m; i++) lmt[i] = jL[i]/info->hxf; 298 } 299 } 300 void Limit3_Fromm(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,const PetscInt sm,const PetscInt mf,const PetscInt fm,const PetscInt ms,PetscInt n,PetscScalar *lmt) 301 { 302 PetscInt i; 303 if (n < sm-1 || n > ms) { 304 for (i=0; i<info->m; i++) lmt[i] = 0.5*(jL[i]+jR[i])/info->hxs; 305 } else if (n == sm-1) { 306 for (i=0; i<info->m; i++) lmt[i] = 0.5*(jL[i]/info->hxs+jR[i]/(info->hxs/2.0+info->hxf/2.0)); 307 } else if (n == sm) { 308 for (i=0; i<info->m; i++) lmt[i] = 0.5*(jL[i]/(info->hxs/2.0+info->hxm/2.0)+jR[i]/info->hxm); 309 } else if (n == ms-1) { 310 for (i=0; i<info->m; i++) lmt[i] = 0.5*(jL[i]/info->hxm+jR[i]/(info->hxs/2.0+info->hxf/2.0)); 311 } else if (n == ms) { 312 for (i=0; i<info->m; i++) lmt[i] = 0.5*(jL[i]/(info->hxm/2.0+info->hxs/2.0)+jR[i]/info->hxs); 313 } else if (n < mf-1 || n > fm) { 314 for (i=0; i<info->m; i++) lmt[i] = 0.5*(jL[i]+jR[i])/info->hxm; 315 } else if (n == mf-1) { 316 for (i=0; i<info->m; i++) lmt[i] = 0.5*(jL[i]/info->hxm+jR[i]/(info->hxm/2.0+info->hxf/2.0)); 317 } else if (n == mf) { 318 for (i=0; i<info->m; i++) lmt[i] = 0.5*(jL[i]/(info->hxm/2.0+info->hxf/2.0)+jR[i]/info->hxf); 319 } else if (n == fm-1) { 320 for (i=0; i<info->m; i++) lmt[i] = 0.5*(jL[i]/info->hxf+jR[i]/(info->hxf/2.0+info->hxm/2.0)); 321 } else if (n == fm) { 322 for (i=0; i<info->m; i++) lmt[i] = 0.5*(jL[i]/(info->hxf/2.0+info->hxm/2.0)+jR[i]/info->hxm); 323 } else { 324 for (i=0; i<info->m; i++) lmt[i] = 0.5*(jL[i]+jR[i])/info->hxf; 325 } 326 } 327 void Limit3_Minmod(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,const PetscInt sm,const PetscInt mf,const PetscInt fm,const PetscInt ms,PetscInt n,PetscScalar *lmt) 328 { 329 PetscInt i; 330 if (n < sm-1 || n > ms) { 331 for (i=0; i<info->m; i++) lmt[i] = MinMod2(jL[i],jR[i])/info->hxs; 332 } else if (n == sm-1) { 333 for (i=0; i<info->m; i++) lmt[i] = MinMod2(jL[i]/info->hxs,jR[i]/(info->hxs/2.0+info->hxf/2.0)); 334 } else if (n == sm) { 335 for (i=0; i<info->m; i++) lmt[i] = MinMod2(jL[i]/(info->hxs/2.0+info->hxf/2.0),jR[i]/info->hxf); 336 } else if (n == ms-1) { 337 for (i=0; i<info->m; i++) lmt[i] = MinMod2(jL[i]/info->hxm,jR[i]/(info->hxm/2.0+info->hxs/2.0)); 338 } else if (n == ms) { 339 for (i=0; i<info->m; i++) lmt[i] = MinMod2(jL[i]/(info->hxm/2.0+info->hxs/2.0),jR[i]/info->hxs); 340 } else if (n < mf-1 || n > fm) { 341 for (i=0; i<info->m; i++) lmt[i] = MinMod2(jL[i],jR[i])/info->hxm; 342 } else if (n == mf-1) { 343 for (i=0; i<info->m; i++) lmt[i] = MinMod2(jL[i]/info->hxm,jR[i]/(info->hxm/2.0+info->hxf/2.0)); 344 } else if (n == mf) { 345 for (i=0; i<info->m; i++) lmt[i] = MinMod2(jL[i]/(info->hxm/2.0+info->hxf/2.0),jR[i]/info->hxf); 346 } else if (n == fm-1) { 347 for (i=0; i<info->m; i++) lmt[i] = MinMod2(jL[i]/info->hxf,jR[i]/(info->hxf/2.0+info->hxm/2.0)); 348 } else if (n == fm) { 349 for (i=0; i<info->m; i++) lmt[i] = MinMod2(jL[i]/(info->hxf/2.0+info->hxm/2.0),jR[i]/info->hxm); 350 } else { 351 for (i=0; i<info->m; i++) lmt[i] = 0.5*(jL[i]+jR[i])/info->hxf; 352 } 353 } 354 void Limit3_Superbee(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,const PetscInt sm,const PetscInt mf,const PetscInt fm,const PetscInt ms,PetscInt n,PetscScalar *lmt) 355 { 356 PetscInt i; 357 if (n < sm-1 || n > ms) { 358 for (i=0; i<info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i],2*jR[i]),MinMod2(2*jL[i],jR[i]))/info->hxs; 359 } else if (n == sm-1) { 360 for (i=0; i<info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i]/info->hxs,2*jR[i]/(info->hxs/2.0+info->hxm/2.0)),MinMod2(2*jL[i]/info->hxs,jR[i]/(info->hxs/2.0+info->hxm/2.0))); 361 } else if (n == sm) { 362 for (i=0; i<info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i]/(info->hxs/2.0+info->hxm/2.0),2*jR[i]/info->hxm),MinMod2(2*jL[i]/(info->hxs/2.0+info->hxm/2.0),jR[i]/info->hxm)); 363 } else if (n == ms-1) { 364 for (i=0; i<info->m; i++) lmt[i] = MinMod2(MinMod2(jL[i]/info->hxm,2*jR[i]/(info->hxm/2.0+info->hxs/2.0)),MinMod2(2*jL[i]/info->hxm,jR[i]/(info->hxm/2.0+info->hxs/2.0))); 365 } else if (n == ms) { 366 for (i=0; i<info->m; i++) lmt[i] = MinMod2(MinMod2(jL[i]/(info->hxm/2.0+info->hxs/2.0),2*jR[i]/info->hxs),MinMod2(2*jL[i]/(info->hxm/2.0+info->hxs/2.0),jR[i]/info->hxs)); 367 } else if (n < mf-1 || n > fm) { 368 for (i=0; i<info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i],2*jR[i]),MinMod2(2*jL[i],jR[i]))/info->hxm; 369 } else if (n == mf-1) { 370 for (i=0; i<info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i]/info->hxm,2*jR[i]/(info->hxm/2.0+info->hxf/2.0)),MinMod2(2*jL[i]/info->hxm,jR[i]/(info->hxm/2.0+info->hxf/2.0))); 371 } else if (n == mf) { 372 for (i=0; i<info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i]/(info->hxm/2.0+info->hxf/2.0),2*jR[i]/info->hxf),MinMod2(2*jL[i]/(info->hxm/2.0+info->hxf/2.0),jR[i]/info->hxf)); 373 } else if (n == fm-1) { 374 for (i=0; i<info->m; i++) lmt[i] = MinMod2(MinMod2(jL[i]/info->hxf,2*jR[i]/(info->hxf/2.0+info->hxm/2.0)),MinMod2(2*jL[i]/info->hxf,jR[i]/(info->hxf/2.0+info->hxm/2.0))); 375 } else if (n == fm) { 376 for (i=0; i<info->m; i++) lmt[i] = MinMod2(MinMod2(jL[i]/(info->hxf/2.0+info->hxm/2.0),2*jR[i]/info->hxm),MinMod2(2*jL[i]/(info->hxf/2.0+info->hxm/2.0),jR[i]/info->hxm)); 377 } else { 378 for (i=0; i<info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i],2*jR[i]),MinMod2(2*jL[i],jR[i]))/info->hxf; 379 } 380 } 381 void Limit3_MC(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,const PetscInt sm,const PetscInt mf,const PetscInt fm,const PetscInt ms,PetscInt n,PetscScalar *lmt) 382 { 383 PetscInt i; 384 if (n < sm-1 || n > ms) { 385 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i],0.5*(jL[i]+jR[i]),2*jR[i])/info->hxs; 386 } else if (n == sm-1) { 387 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i]/info->hxs,0.5*(jL[i]/info->hxs+jR[i]/(info->hxs/2.0+info->hxm/2.0)),2*jR[i]/(info->hxs/2.0+info->hxm/2.0)); 388 } else if (n == sm) { 389 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i]/(info->hxs/2.0+info->hxm/2.0),0.5*(jL[i]/(info->hxs/2.0+info->hxm/2.0)+jR[i]/info->hxm),2*jR[i]/info->hxm); 390 } else if (n == ms-1) { 391 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i]/info->hxm,0.5*(jL[i]/info->hxm+jR[i]/(info->hxm/2.0+info->hxs/2.0)),2*jR[i]/(info->hxm/2.0+info->hxs/2.0)); 392 } else if (n == ms) { 393 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i]/(info->hxm/2.0+info->hxs/2.0),0.5*(jL[i]/(info->hxm/2.0+info->hxs/2.0)+jR[i]/info->hxs),2*jR[i]/info->hxs); 394 } else if (n < mf-1 || n > fm) { 395 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i],0.5*(jL[i]+jR[i]),2*jR[i])/info->hxm; 396 } else if (n == mf-1) { 397 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i]/info->hxm,0.5*(jL[i]/info->hxm+jR[i]/(info->hxm/2.0+info->hxf/2.0)),2*jR[i]/(info->hxm/2.0+info->hxf/2.0)); 398 } else if (n == mf) { 399 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i]/(info->hxm/2.0+info->hxf/2.0),0.5*(jL[i]/(info->hxm/2.0+info->hxf/2.0)+jR[i]/info->hxf),2*jR[i]/info->hxf); 400 } else if (n == fm-1) { 401 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i]/info->hxf,0.5*(jL[i]/info->hxf+jR[i]/(info->hxf/2.0+info->hxm/2.0)),2*jR[i]/(info->hxf/2.0+info->hxm/2.0)); 402 } else if (n == fm) { 403 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i]/(info->hxf/2.0+info->hxm/2.0),0.5*(jL[i]/(info->hxf/2.0+info->hxm/2.0)+jR[i]/info->hxm),2*jR[i]/info->hxm); 404 } else { 405 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i],0.5*(jL[i]+jR[i]),2*jR[i])/info->hxf; 406 } 407 } 408 void Limit3_Koren3(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,const PetscInt sm,const PetscInt mf,const PetscInt fm,const PetscInt ms,PetscInt n,PetscScalar *lmt) 409 { /* Eq 11 of Cada-Torrilhon 2009 */ 410 PetscInt i; 411 if (n < sm-1 || n > ms) { 412 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i],(jL[i]+2*jR[i])/3,2*jR[i])/info->hxs; 413 } else if (n == sm-1) { 414 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i]/info->hxs,(jL[i]/info->hxs+2*jR[i]/(info->hxm/2.0+info->hxs/2.0))/3,2*jR[i]/(info->hxm/2.0+info->hxs/2.0)); 415 } else if (n == sm) { 416 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i]/(info->hxs/2.0+info->hxm/2.0),(jL[i]/(info->hxs/2.0+info->hxm/2.0)+2*jR[i]/info->hxm)/3,2*jR[i]/info->hxm); 417 } else if (n == ms-1) { 418 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i]/info->hxm,(jL[i]/info->hxm+2*jR[i]/(info->hxm/2.0+info->hxs/2.0))/3,2*jR[i]/(info->hxm/2.0+info->hxs/2.0)); 419 } else if (n == ms) { 420 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i]/(info->hxm/2.0+info->hxs/2.0),(jL[i]/(info->hxm/2.0+info->hxs/2.0)+2*jR[i]/info->hxs)/3,2*jR[i]/info->hxs); 421 } else if (n < mf-1 || n > fm) { 422 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i],(jL[i]+2*jR[i])/3,2*jR[i])/info->hxm; 423 } else if (n == mf-1) { 424 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i]/info->hxm,(jL[i]/info->hxm+2*jR[i]/(info->hxm/2.0+info->hxf/2.0))/3,2*jR[i]/(info->hxm/2.0+info->hxf/2.0)); 425 } else if (n == mf) { 426 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i]/(info->hxm/2.0+info->hxf/2.0),(jL[i]/(info->hxm/2.0+info->hxf/2.0)+2*jR[i]/info->hxf)/3,2*jR[i]/info->hxf); 427 } else if (n == fm-1) { 428 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i]/info->hxf,(jL[i]/info->hxf+2*jR[i]/(info->hxf/2.0+info->hxm/2.0))/3,2*jR[i]/(info->hxf/2.0+info->hxm/2.0)); 429 } else if (n == fm) { 430 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i]/(info->hxf/2.0+info->hxm/2.0),(jL[i]/(info->hxf/2.0+info->hxm/2.0)+2*jR[i]/info->hxm)/3,2*jR[i]/info->hxm); 431 } else { 432 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i],(jL[i]+2*jR[i])/3,2*jR[i])/info->hxs; 433 } 434 } 435 PetscErrorCode RiemannListAdd(PetscFunctionList *flist,const char *name,RiemannFunction rsolve) 436 { 437 PetscErrorCode ierr; 438 439 PetscFunctionBeginUser; 440 ierr = PetscFunctionListAdd(flist,name,rsolve);CHKERRQ(ierr); 441 PetscFunctionReturn(0); 442 } 443 444 PetscErrorCode RiemannListFind(PetscFunctionList flist,const char *name,RiemannFunction *rsolve) 445 { 446 PetscErrorCode ierr; 447 448 PetscFunctionBeginUser; 449 ierr = PetscFunctionListFind(flist,name,rsolve);CHKERRQ(ierr); 450 if (!*rsolve) SETERRQ1(PETSC_COMM_SELF,1,"Riemann solver \"%s\" could not be found",name); 451 PetscFunctionReturn(0); 452 } 453 454 PetscErrorCode ReconstructListAdd(PetscFunctionList *flist,const char *name,ReconstructFunction r) 455 { 456 PetscErrorCode ierr; 457 458 PetscFunctionBeginUser; 459 ierr = PetscFunctionListAdd(flist,name,r);CHKERRQ(ierr); 460 PetscFunctionReturn(0); 461 } 462 463 PetscErrorCode ReconstructListFind(PetscFunctionList flist,const char *name,ReconstructFunction *r) 464 { 465 PetscErrorCode ierr; 466 467 PetscFunctionBeginUser; 468 ierr = PetscFunctionListFind(flist,name,r);CHKERRQ(ierr); 469 if (!*r) SETERRQ1(PETSC_COMM_SELF,1,"Reconstruction \"%s\" could not be found",name); 470 PetscFunctionReturn(0); 471 } 472 473 /* --------------------------------- Physics ----------------------------------- */ 474 475 PetscErrorCode PhysicsDestroy_SimpleFree(void *vctx) 476 { 477 PetscErrorCode ierr; 478 479 PetscFunctionBeginUser; 480 ierr = PetscFree(vctx);CHKERRQ(ierr); 481 PetscFunctionReturn(0); 482 } 483 484 /* --------------------------------- Finite Volume Solver ----------------------------------- */ 485 486 PetscErrorCode FVRHSFunction(TS ts,PetscReal time,Vec X,Vec F,void *vctx) 487 { 488 FVCtx *ctx = (FVCtx*)vctx; 489 PetscErrorCode ierr; 490 PetscInt i,j,k,Mx,dof,xs,xm; 491 PetscReal hx,cfl_idt = 0; 492 PetscScalar *x,*f,*slope; 493 Vec Xloc; 494 DM da; 495 496 PetscFunctionBeginUser; 497 ierr = TSGetDM(ts,&da);CHKERRQ(ierr); 498 ierr = DMGetLocalVector(da,&Xloc);CHKERRQ(ierr); /* Xloc contains ghost points */ 499 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 */ 500 hx = (ctx->xmax-ctx->xmin)/Mx; 501 ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr); /* X is solution vector which does not contain ghost points */ 502 ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr); 503 ierr = VecZeroEntries(F);CHKERRQ(ierr); /* F is the right hand side function corresponds to center points */ 504 ierr = DMDAVecGetArray(da,Xloc,&x);CHKERRQ(ierr); 505 ierr = DMDAVecGetArray(da,F,&f);CHKERRQ(ierr); 506 ierr = DMDAGetArray(da,PETSC_TRUE,&slope);CHKERRQ(ierr); /* contains ghost points */ 507 ierr = DMDAGetCorners(da,&xs,0,0,&xm,0,0);CHKERRQ(ierr); 508 509 if (ctx->bctype == FVBC_OUTFLOW) { 510 for (i=xs-2; i<0; i++) { 511 for (j=0; j<dof; j++) x[i*dof+j] = x[j]; 512 } 513 for (i=Mx; i<xs+xm+2; i++) { 514 for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j]; 515 } 516 } 517 for (i=xs-1; i<xs+xm+1; i++) { 518 struct _LimitInfo info; 519 PetscScalar *cjmpL,*cjmpR; 520 /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */ 521 ierr = (*ctx->physics.characteristic)(ctx->physics.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds,ctx->xmin+hx*i);CHKERRQ(ierr); 522 /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */ 523 ierr = PetscArrayzero(ctx->cjmpLR,2*dof);CHKERRQ(ierr); 524 cjmpL = &ctx->cjmpLR[0]; 525 cjmpR = &ctx->cjmpLR[dof]; 526 for (j=0; j<dof; j++) { 527 PetscScalar jmpL,jmpR; 528 jmpL = x[(i+0)*dof+j]-x[(i-1)*dof+j]; 529 jmpR = x[(i+1)*dof+j]-x[(i+0)*dof+j]; 530 for (k=0; k<dof; k++) { 531 cjmpL[k] += ctx->Rinv[k+j*dof]*jmpL; 532 cjmpR[k] += ctx->Rinv[k+j*dof]*jmpR; 533 } 534 } 535 /* Apply limiter to the left and right characteristic jumps */ 536 info.m = dof; 537 info.hx = hx; 538 (*ctx->limit)(&info,cjmpL,cjmpR,ctx->cslope); 539 for (j=0; j<dof; j++) ctx->cslope[j] /= hx; /* rescale to a slope */ 540 for (j=0; j<dof; j++) { 541 PetscScalar tmp = 0; 542 for (k=0; k<dof; k++) tmp += ctx->R[j+k*dof]*ctx->cslope[k]; 543 slope[i*dof+j] = tmp; 544 } 545 } 546 547 for (i=xs; i<xs+xm+1; i++) { 548 PetscReal maxspeed; 549 PetscScalar *uL,*uR; 550 uL = &ctx->uLR[0]; 551 uR = &ctx->uLR[dof]; 552 for (j=0; j<dof; j++) { 553 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hx/2; 554 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hx/2; 555 } 556 ierr = (*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed,ctx->xmin+hx*i,ctx->xmin,ctx->xmax);CHKERRQ(ierr); 557 cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hx)); /* Max allowable value of 1/Delta t */ 558 if (i > xs) { 559 for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hx; 560 } 561 if (i < xs+xm) { 562 for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hx; 563 } 564 } 565 ierr = DMDAVecRestoreArray(da,Xloc,&x);CHKERRQ(ierr); 566 ierr = DMDAVecRestoreArray(da,F,&f);CHKERRQ(ierr); 567 ierr = DMDARestoreArray(da,PETSC_TRUE,&slope);CHKERRQ(ierr); 568 ierr = DMRestoreLocalVector(da,&Xloc);CHKERRQ(ierr); 569 ierr = MPI_Allreduce(&cfl_idt,&ctx->cfl_idt,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)da));CHKERRQ(ierr); 570 if (0) { 571 /* We need to a way to inform the TS of a CFL constraint, this is a debugging fragment */ 572 PetscReal dt,tnow; 573 ierr = TSGetTimeStep(ts,&dt);CHKERRQ(ierr); 574 ierr = TSGetTime(ts,&tnow);CHKERRQ(ierr); 575 if (dt > 0.5/ctx->cfl_idt) { 576 if (1) { 577 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); 578 } else SETERRQ2(PETSC_COMM_SELF,1,"Stability constraint exceeded, %g > %g",(double)dt,(double)(ctx->cfl/ctx->cfl_idt)); 579 } 580 } 581 PetscFunctionReturn(0); 582 } 583 584 PetscErrorCode FVSample(FVCtx *ctx,DM da,PetscReal time,Vec U) 585 { 586 PetscErrorCode ierr; 587 PetscScalar *u,*uj; 588 PetscInt i,j,k,dof,xs,xm,Mx; 589 590 PetscFunctionBeginUser; 591 if (!ctx->physics.sample) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Physics has not provided a sampling function"); 592 ierr = DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);CHKERRQ(ierr); 593 ierr = DMDAGetCorners(da,&xs,0,0,&xm,0,0);CHKERRQ(ierr); 594 ierr = DMDAVecGetArray(da,U,&u);CHKERRQ(ierr); 595 ierr = PetscMalloc1(dof,&uj);CHKERRQ(ierr); 596 for (i=xs; i<xs+xm; i++) { 597 const PetscReal h = (ctx->xmax-ctx->xmin)/Mx,xi = ctx->xmin+h/2+i*h; 598 const PetscInt N = 200; 599 /* Integrate over cell i using trapezoid rule with N points. */ 600 for (k=0; k<dof; k++) u[i*dof+k] = 0; 601 for (j=0; j<N+1; j++) { 602 PetscScalar xj = xi+h*(j-N/2)/(PetscReal)N; 603 ierr = (*ctx->physics.sample)(ctx->physics.user,ctx->initial,ctx->bctype,ctx->xmin,ctx->xmax,time,xj,uj);CHKERRQ(ierr); 604 for (k=0; k<dof; k++) u[i*dof+k] += ((j==0 || j==N) ? 0.5 : 1.0)*uj[k]/N; 605 } 606 } 607 ierr = DMDAVecRestoreArray(da,U,&u);CHKERRQ(ierr); 608 ierr = PetscFree(uj);CHKERRQ(ierr); 609 PetscFunctionReturn(0); 610 } 611 612 PetscErrorCode SolutionStatsView(DM da,Vec X,PetscViewer viewer) 613 { 614 PetscErrorCode ierr; 615 PetscReal xmin,xmax; 616 PetscScalar sum,tvsum,tvgsum; 617 const PetscScalar *x; 618 PetscInt imin,imax,Mx,i,j,xs,xm,dof; 619 Vec Xloc; 620 PetscBool iascii; 621 622 PetscFunctionBeginUser; 623 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 624 if (iascii) { 625 /* PETSc lacks a function to compute total variation norm (difficult in multiple dimensions), we do it here */ 626 ierr = DMGetLocalVector(da,&Xloc);CHKERRQ(ierr); 627 ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr); 628 ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr); 629 ierr = DMDAVecGetArrayRead(da,Xloc,(void*)&x);CHKERRQ(ierr); 630 ierr = DMDAGetCorners(da,&xs,0,0,&xm,0,0);CHKERRQ(ierr); 631 ierr = DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);CHKERRQ(ierr); 632 tvsum = 0; 633 for (i=xs; i<xs+xm; i++) { 634 for (j=0; j<dof; j++) tvsum += PetscAbsScalar(x[i*dof+j]-x[(i-1)*dof+j]); 635 } 636 ierr = MPI_Allreduce(&tvsum,&tvgsum,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)da));CHKERRQ(ierr); 637 ierr = DMDAVecRestoreArrayRead(da,Xloc,(void*)&x);CHKERRQ(ierr); 638 ierr = DMRestoreLocalVector(da,&Xloc);CHKERRQ(ierr); 639 ierr = VecMin(X,&imin,&xmin);CHKERRQ(ierr); 640 ierr = VecMax(X,&imax,&xmax);CHKERRQ(ierr); 641 ierr = VecSum(X,&sum);CHKERRQ(ierr); 642 ierr = PetscViewerASCIIPrintf(viewer,"Solution range [%8.5f,%8.5f] with minimum at %D, mean %8.5f, ||x||_TV %8.5f\n",(double)xmin,(double)xmax,imin,(double)(sum/Mx),(double)(tvgsum/Mx));CHKERRQ(ierr); 643 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type not supported"); 644 PetscFunctionReturn(0); 645 } 646