1 2//--------------------------------------------------------------------- 3//--------------------------------------------------------------------- 4#include "header.h" 5 6void exact_rhs() { 7 8//--------------------------------------------------------------------- 9//--------------------------------------------------------------------- 10 11//--------------------------------------------------------------------- 12// compute the right hand side based on exact solution 13//--------------------------------------------------------------------- 14 15 double dtemp[5], xi, eta, zeta, dtpp; 16 int c, m, i, j, k, ip1, im1, jp1, 17 jm1, km1, kp1; 18#define dtemp(m) dtemp[m-1] 19 20 21//--------------------------------------------------------------------- 22// loop over all cells owned by this node 23//--------------------------------------------------------------------- 24 for (c = 1; c <= ncells; c++) { 25 26//--------------------------------------------------------------------- 27// initialize 28//--------------------------------------------------------------------- 29 for (k = 0; k <= cell_size(3,c)-1; k++) { 30 for (j = 0; j <= cell_size(2,c)-1; j++) { 31 for (i = 0; i <= cell_size(1,c)-1; i++) { 32 for (m = 1; m <= 5; m++) { 33 forcing(m,i,j,k,c) = 0.0e0; 34 } 35 } 36 } 37 } 38 39//--------------------------------------------------------------------- 40// xi-direction flux differences 41//--------------------------------------------------------------------- 42 for (k = start(3,c); k <= cell_size(3,c)-end(3,c)-1; k++) { 43 zeta = (double)(k+cell_low(3,c)) * dnzm1; 44 for (j = start(2,c); j <= cell_size(2,c)-end(2,c)-1; j++) { 45 eta = (double)(j+cell_low(2,c)) * dnym1; 46 47 for (i = -2*(1-start(1,c)); i <= cell_size(1,c)+1-2*end(1,c); i++) { 48 xi = (double)(i+cell_low(1,c)) * dnxm1; 49 50 exact_solution(xi, eta, zeta, dtemp); 51 for (m = 1; m <= 5; m++) { 52 ue(i,m) = dtemp(m); 53 } 54 55 dtpp = 1.0e0 / dtemp(1); 56 57 for (m = 2; m <= 5; m++) { 58 buf(i,m) = dtpp * dtemp(m); 59 } 60 61 cuf(i) = buf(i,2) * buf(i,2); 62 buf(i,1) = cuf(i) + buf(i,3) * buf(i,3) + 63 buf(i,4) * buf(i,4) ; 64 q(i) = 0.5e0*(buf(i,2)*ue(i,2) + buf(i,3)*ue(i,3) + 65 buf(i,4)*ue(i,4)); 66 67 } 68 69 for (i = start(1,c); i <= cell_size(1,c)-end(1,c)-1; i++) { 70 im1 = i-1; 71 ip1 = i+1; 72 73 forcing(1,i,j,k,c) = forcing(1,i,j,k,c) - 74 tx2*( ue(ip1,2)-ue(im1,2) )+ 75 dx1tx1*(ue(ip1,1)-2.0e0*ue(i,1)+ue(im1,1)); 76 77 forcing(2,i,j,k,c) = forcing(2,i,j,k,c) - tx2 * ( 78 (ue(ip1,2)*buf(ip1,2)+c2*(ue(ip1,5)-q(ip1)))- 79 (ue(im1,2)*buf(im1,2)+c2*(ue(im1,5)-q(im1))))+ 80 xxcon1*(buf(ip1,2)-2.0e0*buf(i,2)+buf(im1,2))+ 81 dx2tx1*( ue(ip1,2)-2.0e0* ue(i,2)+ue(im1,2)); 82 83 forcing(3,i,j,k,c) = forcing(3,i,j,k,c) - tx2 * ( 84 ue(ip1,3)*buf(ip1,2)-ue(im1,3)*buf(im1,2))+ 85 xxcon2*(buf(ip1,3)-2.0e0*buf(i,3)+buf(im1,3))+ 86 dx3tx1*( ue(ip1,3)-2.0e0*ue(i,3) +ue(im1,3)); 87 88 forcing(4,i,j,k,c) = forcing(4,i,j,k,c) - tx2*( 89 ue(ip1,4)*buf(ip1,2)-ue(im1,4)*buf(im1,2))+ 90 xxcon2*(buf(ip1,4)-2.0e0*buf(i,4)+buf(im1,4))+ 91 dx4tx1*( ue(ip1,4)-2.0e0* ue(i,4)+ ue(im1,4)); 92 93 forcing(5,i,j,k,c) = forcing(5,i,j,k,c) - tx2*( 94 buf(ip1,2)*(c1*ue(ip1,5)-c2*q(ip1))- 95 buf(im1,2)*(c1*ue(im1,5)-c2*q(im1)))+ 96 0.5e0*xxcon3*(buf(ip1,1)-2.0e0*buf(i,1)+ 97 buf(im1,1))+ 98 xxcon4*(cuf(ip1)-2.0e0*cuf(i)+cuf(im1))+ 99 xxcon5*(buf(ip1,5)-2.0e0*buf(i,5)+buf(im1,5))+ 100 dx5tx1*( ue(ip1,5)-2.0e0* ue(i,5)+ ue(im1,5)); 101 } 102 103//--------------------------------------------------------------------- 104// Fourth-order dissipation 105//--------------------------------------------------------------------- 106 if (start(1,c) > 0) { 107 for (m = 1; m <= 5; m++) { 108 i = 1; 109 forcing(m,i,j,k,c) = forcing(m,i,j,k,c) - dssp * 110 (5.0e0*ue(i,m) - 4.0e0*ue(i+1,m) +ue(i+2,m)); 111 i = 2; 112 forcing(m,i,j,k,c) = forcing(m,i,j,k,c) - dssp * 113 (-4.0e0*ue(i-1,m) + 6.0e0*ue(i,m) - 114 4.0e0*ue(i+1,m) + ue(i+2,m)); 115 } 116 } 117 118 for (i = start(1,c)*3; i <= cell_size(1,c)-3*end(1,c)-1; i++) { 119 for (m = 1; m <= 5; m++) { 120 forcing(m,i,j,k,c) = forcing(m,i,j,k,c) - dssp* 121 (ue(i-2,m) - 4.0e0*ue(i-1,m) + 122 6.0e0*ue(i,m) - 4.0e0*ue(i+1,m) + ue(i+2,m)); 123 } 124 } 125 126 if (end(1,c) > 0) { 127 for (m = 1; m <= 5; m++) { 128 i = cell_size(1,c)-3; 129 forcing(m,i,j,k,c) = forcing(m,i,j,k,c) - dssp * 130 (ue(i-2,m) - 4.0e0*ue(i-1,m) + 131 6.0e0*ue(i,m) - 4.0e0*ue(i+1,m)); 132 i = cell_size(1,c)-2; 133 forcing(m,i,j,k,c) = forcing(m,i,j,k,c) - dssp * 134 (ue(i-2,m) - 4.0e0*ue(i-1,m) + 5.0e0*ue(i,m)); 135 } 136 } 137 138 } 139 } 140 141//--------------------------------------------------------------------- 142// eta-direction flux differences 143//--------------------------------------------------------------------- 144 for (k = start(3,c); k <= cell_size(3,c)-end(3,c)-1; k++) { 145 zeta = (double)(k+cell_low(3,c)) * dnzm1; 146 for (i = start(1,c); i <= cell_size(1,c)-end(1,c)-1; i++) { 147 xi = (double)(i+cell_low(1,c)) * dnxm1; 148 149 for (j = -2*(1-start(2,c)); j <= cell_size(2,c)+1-2*end(2,c); j++) { 150 eta = (double)(j+cell_low(2,c)) * dnym1; 151 152 exact_solution(xi, eta, zeta, dtemp); 153 for (m = 1; m <= 5; m++) { 154 ue(j,m) = dtemp(m); 155 } 156 157 dtpp = 1.0e0/dtemp(1); 158 159 for (m = 2; m <= 5; m++) { 160 buf(j,m) = dtpp * dtemp(m); 161 } 162 163 cuf(j) = buf(j,3) * buf(j,3); 164 buf(j,1) = cuf(j) + buf(j,2) * buf(j,2) + 165 buf(j,4) * buf(j,4); 166 q(j) = 0.5e0*(buf(j,2)*ue(j,2) + buf(j,3)*ue(j,3) + 167 buf(j,4)*ue(j,4)); 168 } 169 170 for (j = start(2,c); j <= cell_size(2,c)-end(2,c)-1; j++) { 171 jm1 = j-1; 172 jp1 = j+1; 173 174 forcing(1,i,j,k,c) = forcing(1,i,j,k,c) - 175 ty2*( ue(jp1,3)-ue(jm1,3) )+ 176 dy1ty1*(ue(jp1,1)-2.0e0*ue(j,1)+ue(jm1,1)); 177 178 forcing(2,i,j,k,c) = forcing(2,i,j,k,c) - ty2*( 179 ue(jp1,2)*buf(jp1,3)-ue(jm1,2)*buf(jm1,3))+ 180 yycon2*(buf(jp1,2)-2.0e0*buf(j,2)+buf(jm1,2))+ 181 dy2ty1*( ue(jp1,2)-2.0* ue(j,2)+ ue(jm1,2)); 182 183 forcing(3,i,j,k,c) = forcing(3,i,j,k,c) - ty2*( 184 (ue(jp1,3)*buf(jp1,3)+c2*(ue(jp1,5)-q(jp1)))- 185 (ue(jm1,3)*buf(jm1,3)+c2*(ue(jm1,5)-q(jm1))))+ 186 yycon1*(buf(jp1,3)-2.0e0*buf(j,3)+buf(jm1,3))+ 187 dy3ty1*( ue(jp1,3)-2.0e0*ue(j,3) +ue(jm1,3)); 188 189 forcing(4,i,j,k,c) = forcing(4,i,j,k,c) - ty2*( 190 ue(jp1,4)*buf(jp1,3)-ue(jm1,4)*buf(jm1,3))+ 191 yycon2*(buf(jp1,4)-2.0e0*buf(j,4)+buf(jm1,4))+ 192 dy4ty1*( ue(jp1,4)-2.0e0*ue(j,4)+ ue(jm1,4)); 193 194 forcing(5,i,j,k,c) = forcing(5,i,j,k,c) - ty2*( 195 buf(jp1,3)*(c1*ue(jp1,5)-c2*q(jp1))- 196 buf(jm1,3)*(c1*ue(jm1,5)-c2*q(jm1)))+ 197 0.5e0*yycon3*(buf(jp1,1)-2.0e0*buf(j,1)+ 198 buf(jm1,1))+ 199 yycon4*(cuf(jp1)-2.0e0*cuf(j)+cuf(jm1))+ 200 yycon5*(buf(jp1,5)-2.0e0*buf(j,5)+buf(jm1,5))+ 201 dy5ty1*(ue(jp1,5)-2.0e0*ue(j,5)+ue(jm1,5)); 202 } 203 204//--------------------------------------------------------------------- 205// Fourth-order dissipation 206//--------------------------------------------------------------------- 207 if (start(2,c) > 0) { 208 for (m = 1; m <= 5; m++) { 209 j = 1; 210 forcing(m,i,j,k,c) = forcing(m,i,j,k,c) - dssp * 211 (5.0e0*ue(j,m) - 4.0e0*ue(j+1,m) +ue(j+2,m)); 212 j = 2; 213 forcing(m,i,j,k,c) = forcing(m,i,j,k,c) - dssp * 214 (-4.0e0*ue(j-1,m) + 6.0e0*ue(j,m) - 215 4.0e0*ue(j+1,m) + ue(j+2,m)); 216 } 217 } 218 219 for (j = start(2,c)*3; j <= cell_size(2,c)-3*end(2,c)-1; j++) { 220 for (m = 1; m <= 5; m++) { 221 forcing(m,i,j,k,c) = forcing(m,i,j,k,c) - dssp* 222 (ue(j-2,m) - 4.0e0*ue(j-1,m) + 223 6.0e0*ue(j,m) - 4.0e0*ue(j+1,m) + ue(j+2,m)); 224 } 225 } 226 227 if (end(2,c) > 0) { 228 for (m = 1; m <= 5; m++) { 229 j = cell_size(2,c)-3; 230 forcing(m,i,j,k,c) = forcing(m,i,j,k,c) - dssp * 231 (ue(j-2,m) - 4.0e0*ue(j-1,m) + 232 6.0e0*ue(j,m) - 4.0e0*ue(j+1,m)); 233 j = cell_size(2,c)-2; 234 forcing(m,i,j,k,c) = forcing(m,i,j,k,c) - dssp * 235 (ue(j-2,m) - 4.0e0*ue(j-1,m) + 5.0e0*ue(j,m)); 236 237 } 238 } 239 240 } 241 } 242 243//--------------------------------------------------------------------- 244// zeta-direction flux differences 245//--------------------------------------------------------------------- 246 for (j = start(2,c); j <= cell_size(2,c)-end(2,c)-1; j++) { 247 eta = (double)(j+cell_low(2,c)) * dnym1; 248 for (i = start(1,c); i <= cell_size(1,c)-end(1,c)-1; i++) { 249 xi = (double)(i+cell_low(1,c)) * dnxm1; 250 251 for (k = -2*(1-start(3,c)); k <= cell_size(3,c)+1-2*end(3,c); k++) { 252 zeta = (double)(k+cell_low(3,c)) * dnzm1; 253 254 exact_solution(xi, eta, zeta, dtemp); 255 for (m = 1; m <= 5; m++) { 256 ue(k,m) = dtemp(m); 257 } 258 259 dtpp = 1.0e0/dtemp(1); 260 261 for (m = 2; m <= 5; m++) { 262 buf(k,m) = dtpp * dtemp(m); 263 } 264 265 cuf(k) = buf(k,4) * buf(k,4); 266 buf(k,1) = cuf(k) + buf(k,2) * buf(k,2) + 267 buf(k,3) * buf(k,3); 268 q(k) = 0.5e0*(buf(k,2)*ue(k,2) + buf(k,3)*ue(k,3) + 269 buf(k,4)*ue(k,4)); 270 } 271 272 for (k = start(3,c); k <= cell_size(3,c)-end(3,c)-1; k++) { 273 km1 = k-1; 274 kp1 = k+1; 275 276 forcing(1,i,j,k,c) = forcing(1,i,j,k,c) - 277 tz2*( ue(kp1,4)-ue(km1,4) )+ 278 dz1tz1*(ue(kp1,1)-2.0e0*ue(k,1)+ue(km1,1)); 279 280 forcing(2,i,j,k,c) = forcing(2,i,j,k,c) - tz2 * ( 281 ue(kp1,2)*buf(kp1,4)-ue(km1,2)*buf(km1,4))+ 282 zzcon2*(buf(kp1,2)-2.0e0*buf(k,2)+buf(km1,2))+ 283 dz2tz1*( ue(kp1,2)-2.0e0* ue(k,2)+ ue(km1,2)); 284 285 forcing(3,i,j,k,c) = forcing(3,i,j,k,c) - tz2 * ( 286 ue(kp1,3)*buf(kp1,4)-ue(km1,3)*buf(km1,4))+ 287 zzcon2*(buf(kp1,3)-2.0e0*buf(k,3)+buf(km1,3))+ 288 dz3tz1*(ue(kp1,3)-2.0e0*ue(k,3)+ue(km1,3)); 289 290 forcing(4,i,j,k,c) = forcing(4,i,j,k,c) - tz2 * ( 291 (ue(kp1,4)*buf(kp1,4)+c2*(ue(kp1,5)-q(kp1)))- 292 (ue(km1,4)*buf(km1,4)+c2*(ue(km1,5)-q(km1))))+ 293 zzcon1*(buf(kp1,4)-2.0e0*buf(k,4)+buf(km1,4))+ 294 dz4tz1*( ue(kp1,4)-2.0e0*ue(k,4) +ue(km1,4)); 295 296 forcing(5,i,j,k,c) = forcing(5,i,j,k,c) - tz2 * ( 297 buf(kp1,4)*(c1*ue(kp1,5)-c2*q(kp1))- 298 buf(km1,4)*(c1*ue(km1,5)-c2*q(km1)))+ 299 0.5e0*zzcon3*(buf(kp1,1)-2.0e0*buf(k,1) 300 +buf(km1,1))+ 301 zzcon4*(cuf(kp1)-2.0e0*cuf(k)+cuf(km1))+ 302 zzcon5*(buf(kp1,5)-2.0e0*buf(k,5)+buf(km1,5))+ 303 dz5tz1*( ue(kp1,5)-2.0e0*ue(k,5)+ ue(km1,5)); 304 } 305 306//--------------------------------------------------------------------- 307// Fourth-order dissipation 308//--------------------------------------------------------------------- 309 if (start(3,c) > 0) { 310 for (m = 1; m <= 5; m++) { 311 k = 1; 312 forcing(m,i,j,k,c) = forcing(m,i,j,k,c) - dssp * 313 (5.0e0*ue(k,m) - 4.0e0*ue(k+1,m) +ue(k+2,m)); 314 k = 2; 315 forcing(m,i,j,k,c) = forcing(m,i,j,k,c) - dssp * 316 (-4.0e0*ue(k-1,m) + 6.0e0*ue(k,m) - 317 4.0e0*ue(k+1,m) + ue(k+2,m)); 318 } 319 } 320 321 for (k = start(3,c)*3; k <= cell_size(3,c)-3*end(3,c)-1; k++) { 322 for (m = 1; m <= 5; m++) { 323 forcing(m,i,j,k,c) = forcing(m,i,j,k,c) - dssp* 324 (ue(k-2,m) - 4.0e0*ue(k-1,m) + 325 6.0e0*ue(k,m) - 4.0e0*ue(k+1,m) + ue(k+2,m)); 326 } 327 } 328 329 if (end(3,c) > 0) { 330 for (m = 1; m <= 5; m++) { 331 k = cell_size(3,c)-3; 332 forcing(m,i,j,k,c) = forcing(m,i,j,k,c) - dssp * 333 (ue(k-2,m) - 4.0e0*ue(k-1,m) + 334 6.0e0*ue(k,m) - 4.0e0*ue(k+1,m)); 335 k = cell_size(3,c)-2; 336 forcing(m,i,j,k,c) = forcing(m,i,j,k,c) - dssp * 337 (ue(k-2,m) - 4.0e0*ue(k-1,m) + 5.0e0*ue(k,m)); 338 } 339 } 340 341 } 342 } 343 344//--------------------------------------------------------------------- 345// now change the sign of the forcing function, 346//--------------------------------------------------------------------- 347 for (k = start(3,c); k <= cell_size(3,c)-end(3,c)-1; k++) { 348 for (j = start(2,c); j <= cell_size(2,c)-end(2,c)-1; j++) { 349 for (i = start(1,c); i <= cell_size(1,c)-end(1,c)-1; i++) { 350 for (m = 1; m <= 5; m++) { 351 forcing(m,i,j,k,c) = -1.e0 * forcing(m,i,j,k,c); 352 } 353 } 354 } 355 } 356 357 } 358 359 return; 360} 361