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 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 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 /* ----------------------- 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 PetscErrorCode ierr; 436 437 PetscFunctionBeginUser; 438 ierr = PetscFunctionListAdd(flist,name,rsolve);CHKERRQ(ierr); 439 PetscFunctionReturn(0); 440 } 441 442 PetscErrorCode RiemannListFind(PetscFunctionList flist,const char *name,RiemannFunction *rsolve) 443 { 444 PetscErrorCode ierr; 445 446 PetscFunctionBeginUser; 447 ierr = PetscFunctionListFind(flist,name,rsolve);CHKERRQ(ierr); 448 if (!*rsolve) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Riemann solver \"%s\" could not be found",name); 449 PetscFunctionReturn(0); 450 } 451 452 PetscErrorCode ReconstructListAdd(PetscFunctionList *flist,const char *name,ReconstructFunction r) 453 { 454 PetscErrorCode ierr; 455 456 PetscFunctionBeginUser; 457 ierr = PetscFunctionListAdd(flist,name,r);CHKERRQ(ierr); 458 PetscFunctionReturn(0); 459 } 460 461 PetscErrorCode ReconstructListFind(PetscFunctionList flist,const char *name,ReconstructFunction *r) 462 { 463 PetscErrorCode ierr; 464 465 PetscFunctionBeginUser; 466 ierr = PetscFunctionListFind(flist,name,r);CHKERRQ(ierr); 467 if (!*r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Reconstruction \"%s\" could not be found",name); 468 PetscFunctionReturn(0); 469 } 470 471 PetscErrorCode RiemannListAdd_2WaySplit(PetscFunctionList *flist,const char *name,RiemannFunction_2WaySplit rsolve) 472 { 473 PetscErrorCode ierr; 474 475 PetscFunctionBeginUser; 476 ierr = PetscFunctionListAdd(flist,name,rsolve);CHKERRQ(ierr); 477 PetscFunctionReturn(0); 478 } 479 480 PetscErrorCode RiemannListFind_2WaySplit(PetscFunctionList flist,const char *name,RiemannFunction_2WaySplit *rsolve) 481 { 482 PetscErrorCode ierr; 483 484 PetscFunctionBeginUser; 485 ierr = PetscFunctionListFind(flist,name,rsolve);CHKERRQ(ierr); 486 if (!*rsolve) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Riemann solver \"%s\" could not be found",name); 487 PetscFunctionReturn(0); 488 } 489 490 PetscErrorCode ReconstructListAdd_2WaySplit(PetscFunctionList *flist,const char *name,ReconstructFunction_2WaySplit r) 491 { 492 PetscErrorCode ierr; 493 494 PetscFunctionBeginUser; 495 ierr = PetscFunctionListAdd(flist,name,r);CHKERRQ(ierr); 496 PetscFunctionReturn(0); 497 } 498 499 PetscErrorCode ReconstructListFind_2WaySplit(PetscFunctionList flist,const char *name,ReconstructFunction_2WaySplit *r) 500 { 501 PetscErrorCode ierr; 502 503 PetscFunctionBeginUser; 504 ierr = PetscFunctionListFind(flist,name,r);CHKERRQ(ierr); 505 if (!*r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Reconstruction \"%s\" could not be found",name); 506 PetscFunctionReturn(0); 507 } 508 509 /* --------------------------------- Physics ------- */ 510 PetscErrorCode PhysicsDestroy_SimpleFree(void *vctx) 511 { 512 PetscErrorCode ierr; 513 514 PetscFunctionBeginUser; 515 ierr = PetscFree(vctx);CHKERRQ(ierr); 516 PetscFunctionReturn(0); 517 } 518 519 /* --------------------------------- Finite Volume Solver --------------- */ 520 PetscErrorCode FVRHSFunction(TS ts,PetscReal time,Vec X,Vec F,void *vctx) 521 { 522 FVCtx *ctx = (FVCtx*)vctx; 523 PetscErrorCode ierr; 524 PetscInt i,j,k,Mx,dof,xs,xm; 525 PetscReal hx,cfl_idt = 0; 526 PetscScalar *x,*f,*slope; 527 Vec Xloc; 528 DM da; 529 530 PetscFunctionBeginUser; 531 ierr = TSGetDM(ts,&da);CHKERRQ(ierr); 532 ierr = DMGetLocalVector(da,&Xloc);CHKERRQ(ierr); /* Xloc contains ghost points */ 533 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 */ 534 hx = (ctx->xmax-ctx->xmin)/Mx; 535 ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr); /* X is solution vector which does not contain ghost points */ 536 ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr); 537 ierr = VecZeroEntries(F);CHKERRQ(ierr); /* F is the right hand side function corresponds to center points */ 538 ierr = DMDAVecGetArray(da,Xloc,&x);CHKERRQ(ierr); 539 ierr = DMDAVecGetArray(da,F,&f);CHKERRQ(ierr); 540 ierr = DMDAGetArray(da,PETSC_TRUE,&slope);CHKERRQ(ierr); /* contains ghost points */ 541 ierr = DMDAGetCorners(da,&xs,0,0,&xm,0,0);CHKERRQ(ierr); 542 543 if (ctx->bctype == FVBC_OUTFLOW) { 544 for (i=xs-2; i<0; i++) { 545 for (j=0; j<dof; j++) x[i*dof+j] = x[j]; 546 } 547 for (i=Mx; i<xs+xm+2; i++) { 548 for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j]; 549 } 550 } 551 552 for (i=xs-1; i<xs+xm+1; i++) { 553 struct _LimitInfo info; 554 PetscScalar *cjmpL,*cjmpR; 555 /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */ 556 ierr = (*ctx->physics.characteristic)(ctx->physics.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds,ctx->xmin+hx*i);CHKERRQ(ierr); 557 /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */ 558 ierr = PetscArrayzero(ctx->cjmpLR,2*dof);CHKERRQ(ierr); 559 cjmpL = &ctx->cjmpLR[0]; 560 cjmpR = &ctx->cjmpLR[dof]; 561 for (j=0; j<dof; j++) { 562 PetscScalar jmpL,jmpR; 563 jmpL = x[(i+0)*dof+j]-x[(i-1)*dof+j]; 564 jmpR = x[(i+1)*dof+j]-x[(i+0)*dof+j]; 565 for (k=0; k<dof; k++) { 566 cjmpL[k] += ctx->Rinv[k+j*dof]*jmpL; 567 cjmpR[k] += ctx->Rinv[k+j*dof]*jmpR; 568 } 569 } 570 /* Apply limiter to the left and right characteristic jumps */ 571 info.m = dof; 572 info.hx = hx; 573 (*ctx->limit)(&info,cjmpL,cjmpR,ctx->cslope); 574 for (j=0; j<dof; j++) ctx->cslope[j] /= hx; /* rescale to a slope */ 575 for (j=0; j<dof; j++) { 576 PetscScalar tmp = 0; 577 for (k=0; k<dof; k++) tmp += ctx->R[j+k*dof]*ctx->cslope[k]; 578 slope[i*dof+j] = tmp; 579 } 580 } 581 582 for (i=xs; i<xs+xm+1; i++) { 583 PetscReal maxspeed; 584 PetscScalar *uL,*uR; 585 uL = &ctx->uLR[0]; 586 uR = &ctx->uLR[dof]; 587 for (j=0; j<dof; j++) { 588 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hx/2; 589 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hx/2; 590 } 591 ierr = (*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed,ctx->xmin+hx*i,ctx->xmin,ctx->xmax);CHKERRQ(ierr); 592 cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hx)); /* Max allowable value of 1/Delta t */ 593 if (i > xs) { 594 for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hx; 595 } 596 if (i < xs+xm) { 597 for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hx; 598 } 599 } 600 ierr = DMDAVecRestoreArray(da,Xloc,&x);CHKERRQ(ierr); 601 ierr = DMDAVecRestoreArray(da,F,&f);CHKERRQ(ierr); 602 ierr = DMDARestoreArray(da,PETSC_TRUE,&slope);CHKERRQ(ierr); 603 ierr = DMRestoreLocalVector(da,&Xloc);CHKERRQ(ierr); 604 ierr = MPI_Allreduce(&cfl_idt,&ctx->cfl_idt,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)da));CHKERRMPI(ierr); 605 if (0) { 606 /* We need to a way to inform the TS of a CFL constraint, this is a debugging fragment */ 607 PetscReal dt,tnow; 608 ierr = TSGetTimeStep(ts,&dt);CHKERRQ(ierr); 609 ierr = TSGetTime(ts,&tnow);CHKERRQ(ierr); 610 if (dt > 0.5/ctx->cfl_idt) { 611 if (1) { 612 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); 613 } else SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Stability constraint exceeded, %g > %g",(double)dt,(double)(ctx->cfl/ctx->cfl_idt)); 614 } 615 } 616 PetscFunctionReturn(0); 617 } 618 619 PetscErrorCode FVSample(FVCtx *ctx,DM da,PetscReal time,Vec U) 620 { 621 PetscErrorCode ierr; 622 PetscScalar *u,*uj; 623 PetscInt i,j,k,dof,xs,xm,Mx; 624 625 PetscFunctionBeginUser; 626 if (!ctx->physics.sample) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Physics has not provided a sampling function"); 627 ierr = DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);CHKERRQ(ierr); 628 ierr = DMDAGetCorners(da,&xs,0,0,&xm,0,0);CHKERRQ(ierr); 629 ierr = DMDAVecGetArray(da,U,&u);CHKERRQ(ierr); 630 ierr = PetscMalloc1(dof,&uj);CHKERRQ(ierr); 631 for (i=xs; i<xs+xm; i++) { 632 const PetscReal h = (ctx->xmax-ctx->xmin)/Mx,xi = ctx->xmin+h/2+i*h; 633 const PetscInt N = 200; 634 /* Integrate over cell i using trapezoid rule with N points. */ 635 for (k=0; k<dof; k++) u[i*dof+k] = 0; 636 for (j=0; j<N+1; j++) { 637 PetscScalar xj = xi+h*(j-N/2)/(PetscReal)N; 638 ierr = (*ctx->physics.sample)(ctx->physics.user,ctx->initial,ctx->bctype,ctx->xmin,ctx->xmax,time,xj,uj);CHKERRQ(ierr); 639 for (k=0; k<dof; k++) u[i*dof+k] += ((j==0 || j==N) ? 0.5 : 1.0)*uj[k]/N; 640 } 641 } 642 ierr = DMDAVecRestoreArray(da,U,&u);CHKERRQ(ierr); 643 ierr = PetscFree(uj);CHKERRQ(ierr); 644 PetscFunctionReturn(0); 645 } 646 647 PetscErrorCode SolutionStatsView(DM da,Vec X,PetscViewer viewer) 648 { 649 PetscErrorCode ierr; 650 PetscReal xmin,xmax; 651 PetscScalar sum,tvsum,tvgsum; 652 const PetscScalar *x; 653 PetscInt imin,imax,Mx,i,j,xs,xm,dof; 654 Vec Xloc; 655 PetscBool iascii; 656 657 PetscFunctionBeginUser; 658 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 659 if (iascii) { 660 /* PETSc lacks a function to compute total variation norm (difficult in multiple dimensions), we do it here */ 661 ierr = DMGetLocalVector(da,&Xloc);CHKERRQ(ierr); 662 ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr); 663 ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr); 664 ierr = DMDAVecGetArrayRead(da,Xloc,(void*)&x);CHKERRQ(ierr); 665 ierr = DMDAGetCorners(da,&xs,0,0,&xm,0,0);CHKERRQ(ierr); 666 ierr = DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);CHKERRQ(ierr); 667 tvsum = 0; 668 for (i=xs; i<xs+xm; i++) { 669 for (j=0; j<dof; j++) tvsum += PetscAbsScalar(x[i*dof+j]-x[(i-1)*dof+j]); 670 } 671 ierr = MPI_Allreduce(&tvsum,&tvgsum,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)da));CHKERRMPI(ierr); 672 ierr = DMDAVecRestoreArrayRead(da,Xloc,(void*)&x);CHKERRQ(ierr); 673 ierr = DMRestoreLocalVector(da,&Xloc);CHKERRQ(ierr); 674 ierr = VecMin(X,&imin,&xmin);CHKERRQ(ierr); 675 ierr = VecMax(X,&imax,&xmax);CHKERRQ(ierr); 676 ierr = VecSum(X,&sum);CHKERRQ(ierr); 677 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); 678 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type not supported"); 679 PetscFunctionReturn(0); 680 } 681