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