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