最近接到一个任务，就是解决FVCOM中对流扩散计算不守衡问题。导师认为是其求解时候水平和垂向计算分开求解所导致的，目前我也没搞清到底有什么问题，反正就是让把SUNTANS的对流扩散计算挪到FVCOM中。
SUNTANS模型手册：http://web.stanford.edu/group/suntans/cgi-bin/documentation/user_guide/user_guide.html
介绍文献：《An unstructured-grid, finite-volume, nonhydrostatic, parallel coastal ocean simulator》
代码所谓研究讨论之用这里只公布部分：
1 /* 2 * File: scalars.c 3 * Author: Oliver B. Fringer 4 * Institution: Stanford University 5 * ---------------------------------------- 6 * This file contains the scalar transport function. 7 * 8 * Copyright (C) 2005-2006 The Board of Trustees of the Leland Stanford Junior 9 * University. All Rights Reserved. 10 * 11 */ 12 #include "scalars.h" 13 #include "util.h" 14 #include "tvd.h" 15 #include "initialization.h" 16 17 #define SMALL_CONSISTENCY 1e-5 18 19 REAL smin_value, smax_value; 20 21 /* 22 * Function: UpdateScalars 23 * Usage: UpdateScalars(grid,phys,prop,wnew,scalar,Cn,kappa,kappaH,kappa_tv,theta); 24 * --------------------------------------------------------------------------- 25 * Update the scalar quantity stored in the array denoted by scal using the 26 * theta method for vertical advection and vertical diffusion and Adams-Bashforth 27 * for horizontal advection and diffusion. 28 * 29 * Cn must store the AB terms from time step n-1 for this scalar 30 * kappa denotes the vertical scalar diffusivity 31 * kappaH denotes the horizontal scalar diffusivity 32 * kappa_tv denotes the vertical turbulent scalar diffusivity 33 * 34 */ 35 void UpdateScalars(gridT *grid, physT *phys, propT *prop, REAL **wnew, REAL **scal, REAL **boundary_scal, REAL **Cn, 36 REAL kappa, REAL kappaH, REAL **kappa_tv, REAL theta, 37 REAL **src1, REAL **src2, REAL *Ftop, REAL *Fbot, int alpha_top, int alpha_bot, 38 MPI_Comm comm, int myproc, int checkflag, int TVDscheme) 39 { 40 int i, iptr, j, jptr, ib, k, nf, ktop; 41 int Nc=grid->Nc, normal, nc1, nc2, ne; 42 REAL df, dg, Ac, dt=prop->dt, fab, *a, *b, *c, *d, *ap, *am, *bd, *uflux, dznew, mass, *sp, *temp; 43 REAL smin, smax, div_local, div_da; 44 int k1, k2, kmin, imin, kmax, imax, mincount, maxcount, allmincount, allmaxcount, flag; 45 46 prop->TVD = TVDscheme; 47 // These are used mostly debugging to turn on/off vertical and horizontal TVD. 48 prop->horiTVD = 1; 49 prop->vertTVD = 1; 50 51 ap = phys->ap; 52 am = phys->am; 53 bd = phys->bp; 54 temp = phys->bm; 55 a = phys->a; 56 b = phys->b; 57 c = phys->c; 58 d = phys->d; 59 60 // Never use AB2 61 if(1) { 62 fab=1; 63 for(i=0;i<grid->Nc;i++) 64 for(k=0;k<grid->Nk[i];k++) 65 Cn[i][k]=0; 66 } else 67 fab=1.5; 68 69 for(i=0;i<Nc;i++) 70 for(k=0;k<grid->Nk[i];k++) 71 phys->stmp[i][k]=scal[i][k]; 72 73 // Add on boundary fluxes, using stmp2 as the temporary storage 74 // variable 75 //for(iptr=grid->celldist[0];iptr<grid->celldist[1];iptr++) { 76 for(iptr=grid->celldist[0];iptr<grid->celldist[2];iptr++) { 77 i = grid->cellp[iptr]; 78 79 for(k=grid->ctop[i];k<grid->Nk[i];k++) 80 phys->stmp2[i][k]=0; 81 } 82 83 if(boundary_scal) { 84 for(jptr=grid->edgedist[2];jptr<grid->edgedist[5];jptr++) { 85 j = grid->edgep[jptr]; 86 87 ib = grid->grad[2*j]; 88 89 // Set the value of stmp2 adjacent to the boundary to the value of the boundary. 90 // This will be used to add the boundary flux when stmp2 is used again below. 91 for(k=grid->ctop[ib];k<grid->Nk[ib];k++) 92 phys->stmp2[ib][k]=boundary_scal[jptr-grid->edgedist[2]][k]; 93 } 94 } 95 96 // Compute the scalar on the vertical faces (for horiz. advection) 97 98 if(prop->TVD && prop->horiTVD) 99 HorizontalFaceScalars(grid,phys,prop,scal,boundary_scal,prop->TVD,comm,myproc); 100 101 //for(iptr=grid->celldist[0];iptr<grid->celldist[1];iptr++) { 102 for(iptr=grid->celldist[0];iptr<grid->celldist[2];iptr++) { 103 i = grid->cellp[iptr]; 104 Ac = grid->Ac[i]; 105 106 if(grid->ctop[i]>=grid->ctopold[i]) { 107 ktop=grid->ctop[i]; 108 dznew=grid->dzz[i][ktop]; 109 } else { 110 ktop=grid->ctopold[i]; 111 dznew=0; 112 for(k=grid->ctop[i];k<=grid->ctopold[i];k++) 113 dznew+=grid->dzz[i][k]; 114 } 115 116 // These are the advective components of the tridiagonal 117 // at the new time step. 118 if(!(prop->TVD && prop->vertTVD)) 119 for(k=0;k<grid->Nk[i]+1;k++) { 120 ap[k] = 0.5*(wnew[i][k]+fabs(wnew[i][k])); 121 am[k] = 0.5*(wnew[i][k]-fabs(wnew[i][k])); 122 } 123 else // Compute the ap/am for TVD schemes 124 GetApAm(ap,am,phys->wp,phys->wm,phys->Cp,phys->Cm,phys->rp,phys->rm, 125 wnew,grid->dzz,scal,i,grid->Nk[i],ktop,prop->dt,prop->TVD); 126 127 for(k=ktop+1;k<grid->Nk[i];k++){ 128 a[k-ktop]=theta*dt*am[k]; 129 b[k-ktop]=grid->dzz[i][k]+theta*dt*(ap[k]-am[k+1]); 130 c[k-ktop]=-theta*dt*ap[k+1]; 131 } 132 133 // Top cell advection 134 a[0]=0; 135 b[0]=dznew-theta*dt*am[ktop+1]; 136 c[0]=-theta*dt*ap[ktop+1]; 137 138 // Bottom cell no-flux boundary condition for advection 139 b[(grid->Nk[i]-1)-ktop]+=c[(grid->Nk[i]-1)-ktop]; 140 141 // Implicit vertical diffusion terms 142 for(k=ktop+1;k<grid->Nk[i];k++) 143 bd[k]=(2.0*kappa+kappa_tv[i][k-1]+kappa_tv[i][k])/ 144 (grid->dzz[i][k-1]+grid->dzz[i][k]); 145 146 for(k=ktop+1;k<grid->Nk[i]-1;k++) { 147 a[k-ktop]-=theta*dt*bd[k]; 148 b[k-ktop]+=theta*dt*(bd[k]+bd[k+1]); 149 c[k-ktop]-=theta*dt*bd[k+1]; 150 } 151 if(src1) 152 for(k=ktop;k<grid->Nk[i];k++) 153 b[k-ktop]+=grid->dzz[i][k]*src1[i][k]*theta*dt; 154 155 // Diffusive fluxes only when more than 1 layer 156 if(ktop<grid->Nk[i]-1) { 157 // Top cell diffusion 158 b[0]+=theta*dt*(bd[ktop+1]+2*alpha_top*bd[ktop+1]); 159 c[0]-=theta*dt*bd[ktop+1]; 160 161 // Bottom cell diffusion 162 a[(grid->Nk[i]-1)-ktop]-=theta*dt*bd[grid->Nk[i]-1]; 163 b[(grid->Nk[i]-1)-ktop]+=theta*dt*(bd[grid->Nk[i]-1]+2*alpha_bot*bd[grid->Nk[i]-1]); 164 } 165 166 // Explicit part into source term d[] 167 for(k=ktop+1;k<grid->Nk[i];k++) 168 d[k-ktop]=grid->dzzold[i][k]*phys->stmp[i][k]; 169 if(src1) 170 for(k=ktop+1;k<grid->Nk[i];k++) 171 d[k-ktop]-=src1[i][k]*(1-theta)*dt*grid->dzzold[i][k]*phys->stmp[i][k]; 172 173 d[0]=0; 174 if(grid->ctopold[i]<=grid->ctop[i]) { 175 for(k=grid->ctopold[i];k<=grid->ctop[i];k++) 176 d[0]+=grid->dzzold[i][k]*phys->stmp[i][k]; 177 if(src1) 178 for(k=grid->ctopold[i];k<=grid->ctop[i];k++) 179 d[0]-=src1[i][k]*(1-theta)*dt*grid->dzzold[i][k]*phys->stmp[i][k]; 180 } else { 181 d[0]=grid->dzzold[i][ktop]*phys->stmp[i][ktop]; 182 if(src1) 183 d[0]-=src1[i][ktop]*(1-theta)*dt*grid->dzzold[i][ktop]*phys->stmp[i][k]; 184 } 185 186 // These are the advective components of the tridiagonal 187 // that use the new velocity 188 if(!(prop->TVD && prop->vertTVD)) 189 for(k=0;k<grid->Nk[i]+1;k++) { 190 ap[k] = 0.5*(phys->wtmp2[i][k]+fabs(phys->wtmp2[i][k])); 191 am[k] = 0.5*(phys->wtmp2[i][k]-fabs(phys->wtmp2[i][k])); 192 } 193 else // Compute the ap/am for TVD schemes 194 GetApAm(ap,am,phys->wp,phys->wm,phys->Cp,phys->Cm,phys->rp,phys->rm, 195 phys->wtmp2,grid->dzzold,phys->stmp,i,grid->Nk[i],ktop,prop->dt,prop->TVD); 196 197 // Explicit advection and diffusion 198 for(k=ktop+1;k<grid->Nk[i]-1;k++) 199 d[k-ktop]-=(1-theta)*dt*(am[k]*phys->stmp[i][k-1]+ 200 (ap[k]-am[k+1])*phys->stmp[i][k]- 201 ap[k+1]*phys->stmp[i][k+1])- 202 (1-theta)*dt*(bd[k]*phys->stmp[i][k-1] 203 -(bd[k]+bd[k+1])*phys->stmp[i][k] 204 +bd[k+1]*phys->stmp[i][k+1]); 205 206 if(ktop<grid->Nk[i]-1) { 207 //Flux through bottom of top cell 208 k=ktop; 209 d[0]=d[0]-(1-theta)*dt*(-am[k+1]*phys->stmp[i][k]- 210 ap[k+1]*phys->stmp[i][k+1])+ 211 (1-theta)*dt*(-(2*alpha_top*bd[k+1]+bd[k+1])*phys->stmp[i][k]+ 212 bd[k+1]*phys->stmp[i][k+1]); 213 if(Ftop) d[0]+=dt*(1-alpha_top+2*alpha_top*bd[k+1])*Ftop[i]; 214 215 // Through top of bottom cell 216 k=grid->Nk[i]-1; 217 d[k-ktop]-=(1-theta)*dt*(am[k]*phys->stmp[i][k-1]+ 218 ap[k]*phys->stmp[i][k])- 219 (1-theta)*dt*(bd[k]*phys->stmp[i][k-1]- 220 (bd[k]+2*alpha_bot*bd[k])*phys->stmp[i][k]); 221 if(Fbot) d[k-ktop]+=dt*(-1+alpha_bot+2*alpha_bot*bd[k])*Fbot[i]; 222 } 223 224 // First add on the source term from the previous time step. 225 if(grid->ctop[i]<=grid->ctopold[i]) { 226 for(k=grid->ctop[i];k<=grid->ctopold[i];k++) 227 d[0]+=(1-fab)*Cn[i][grid->ctopold[i]]/(1+fabs(grid->ctop[i]-grid->ctopold[i])); 228 for(k=grid->ctopold[i]+1;k<grid->Nk[i];k++) 229 d[k-grid->ctopold[i]]+=(1-fab)*Cn[i][k]; 230 } else { 231 for(k=grid->ctopold[i];k<=grid->ctop[i];k++) 232 d[0]+=(1-fab)*Cn[i][k]; 233 for(k=grid->ctop[i]+1;k<grid->Nk[i];k++) 234 d[k-grid->ctop[i]]+=(1-fab)*Cn[i][k]; 235 } 236 237 for(k=0;k<grid->ctop[i];k++) 238 Cn[i][k]=0; 239 240 if(src2) 241 for(k=grid->ctop[i];k<grid->Nk[i];k++) 242 Cn[i][k-ktop]=dt*src2[i][k]*grid->dzzold[i][k]; 243 else 244 for(k=grid->ctop[i];k<grid->Nk[i];k++) 245 Cn[i][k]=0; 246 247 // Now create the source term for the current time step 248 for(k=0;k<grid->Nk[i];k++) 249 ap[k]=0; 250 251 for(nf=0;nf<grid->nfaces[i];nf++) { 252 ne = grid->face[i*grid->maxfaces+nf]; 253 normal = grid->normal[i*grid->maxfaces+nf]; 254 df = grid->df[ne]; 255 dg = grid->dg[ne]; 256 nc1 = grid->grad[2*ne]; 257 nc2 = grid->grad[2*ne+1]; 258 if(nc1==-1) nc1=nc2; 259 if(nc2==-1) { 260 nc2=nc1; 261 if(boundary_scal && (grid->mark[ne]==2 || grid->mark[ne]==3)) 262 sp=phys->stmp2[nc1]; 263 else 264 sp=phys->stmp[nc1]; 265 } else 266 sp=phys->stmp[nc2]; 267 268 if(!(prop->TVD && prop->horiTVD)) { 269 for(k=0;k<grid->Nke[ne];k++) 270 temp[k]=UpWind(phys->utmp2[ne][k], 271 phys->stmp[nc1][k], 272 sp[k]); 273 } else { 274 for(k=0;k<grid->Nke[ne];k++) 275 if(phys->utmp2[ne][k]>0) 276 temp[k]=phys->SfHp[ne][k]; 277 else 278 temp[k]=phys->SfHm[ne][k]; 279 } 280 281 for(k=0;k<grid->Nke[ne];k++) 282 ap[k] += dt*df*normal/Ac*(theta*phys->u[ne][k]+(1-theta)*phys->utmp2[ne][k]) 283 *temp[k]*grid->dzf[ne][k]; 284 } 285 286 for(k=ktop+1;k<grid->Nk[i];k++) 287 Cn[i][k-ktop]-=ap[k]; 288 289 for(k=0;k<=ktop;k++) 290 Cn[i][0]-=ap[k]; 291 292 // Add on the source from the current time step to the rhs. 293 for(k=0;k<grid->Nk[i]-ktop;k++) 294 d[k]+=fab*Cn[i][k]; 295 296 // Add on the volume correction if h was < -d 297 /* 298 if(grid->ctop[i]==grid->Nk[i]-1) 299 d[grid->Nk[i]-ktop-1]+=phys->hcorr[i]*phys->stmp[i][grid->ctop[i]]; 300 */ 301 302 for(k=ktop;k<grid->Nk[i];k++) 303 ap[k]=Cn[i][k-ktop]; 304 for(k=0;k<=ktop;k++) 305 Cn[i][k]=0; 306 for(k=ktop+1;k<grid->Nk[i];k++) 307 Cn[i][k]=ap[k]; 308 for(k=grid->ctop[i];k<=ktop;k++) 309 Cn[i][k]=ap[ktop]/(1+fabs(grid->ctop[i]-ktop)); 310 311 if(grid->Nk[i]-ktop>1) 312 TriSolve(a,b,c,d,&(scal[i][ktop]),grid->Nk[i]-ktop); 313 else if(prop->n>1) { 314 if(b[0]>0 && phys->active[i]) 315 scal[i][ktop]=d[0]/b[0]; 316 else 317 scal[i][ktop]=0; 318 } 319 320 for(k=0;k<grid->ctop[i];k++) 321 scal[i][k]=0; 322 323 for(k=grid->ctop[i];k<grid->ctopold[i];k++) 324 scal[i][k]=scal[i][ktop]; 325 } 326 327 // Code to check divergence change CHECKCONSISTENCY to 1 in suntans.h 328 if(CHECKCONSISTENCY && checkflag) { 329 330 if(prop->n==1+prop->nstart) { 331 smin=INFTY; 332 smax=-INFTY; 333 for(i=0;i<grid->Nc;i++) { 334 for(k=grid->ctop[i];k<grid->Nk[i];k++) { 335 if(phys->stmp[i][k]>smax) { 336 smax=phys->stmp[i][k]; 337 imax=i; 338 kmax=k; 339 } 340 if(phys->stmp[i][k]<smin) { 341 smin=phys->stmp[i][k]; 342 imin=i; 343 kmin=k; 344 } 345 } 346 } 347 MPI_Reduce(&smin,&smin_value,1,MPI_DOUBLE,MPI_MIN,0,comm); 348 MPI_Reduce(&smax,&smax_value,1,MPI_DOUBLE,MPI_MAX,0,comm); 349 MPI_Bcast(&smin_value,1,MPI_DOUBLE,0,comm); 350 MPI_Bcast(&smax_value,1,MPI_DOUBLE,0,comm); 351 352 if(myproc==0) 353 printf("Minimum scalar: %.2f, maximum: %.2f\n",smin_value,smax_value); 354 } 355 356 //for(iptr=grid->celldist[0];iptr<grid->celldist[1];iptr++) { 357 for(iptr=grid->celldist[0];iptr<grid->celldist[2];iptr++) { 358 i = grid->cellp[iptr]; 359 360 flag=0; 361 for(nf=0;nf<grid->nfaces[i];nf++) { 362 if(grid->mark[grid->face[i*grid->maxfaces+nf]]==2 || 363 grid->mark[grid->face[i*grid->maxfaces+nf]]==3) { 364 flag=1; 365 break; 366 } 367 } 368 369 if(!flag) { 370 div_da=0; 371 372 for(k=0;k<grid->Nk[i];k++) { 373 div_da+=grid->Ac[i]*(grid->dzz[i][k]-grid->dzzold[i][k])/prop->dt; 374 375 div_local=0; 376 for(nf=0;nf<grid->nfaces[i];nf++) { 377 ne=grid->face[i*grid->maxfaces+nf]; 378 div_local+=(theta*phys->u[ne][k]+(1-theta)*phys->utmp2[ne][k]) 379 *grid->dzf[ne][k]*grid->normal[i*grid->maxfaces+nf]*grid->df[ne]; 380 } 381 div_da+=div_local; 382 div_local+=grid->Ac[i]*(theta*(wnew[i][k]-wnew[i][k+1])+ 383 (1-theta)*(phys->wtmp2[i][k]-phys->wtmp2[i][k+1])); 384 385 if(k>=grid->ctop[i]) { 386 if(fabs(div_local)>SMALL_CONSISTENCY && grid->dzz[imin][0]>DRYCELLHEIGHT) 387 printf("Step: %d, proc: %d, locally-divergent at %d, %d, div=%e\n", 388 prop->n,myproc,i,k,div_local); 389 } 390 } 391 if(fabs(div_da)>SMALL_CONSISTENCY && phys->h[i]+grid->dv[i]>DRYCELLHEIGHT) 392 printf("Step: %d, proc: %d, Depth-Ave divergent at i=%d, div=%e\n", 393 prop->n,myproc,i,div_da); 394 } 395 } 396 397 mincount=0; 398 maxcount=0; 399 smin=INFTY; 400 smax=-INFTY; 401 //for(iptr=grid->celldist[0];iptr<grid->celldist[1];iptr++) { 402 for(iptr=grid->celldist[0];iptr<grid->celldist[2];iptr++) { 403 i = grid->cellp[iptr]; 404 405 flag=0; 406 for(nf=0;nf<grid->nfaces[i];nf++) { 407 if(grid->mark[grid->face[i*grid->maxfaces+nf]]==2 || grid->mark[grid->face[i*grid->maxfaces+nf]]==3) { 408 flag=1; 409 break; 410 } 411 } 412 413 if(!flag) { 414 for(k=grid->ctop[i];k<grid->Nk[i];k++) { 415 if(scal[i][k]>smax) { 416 smax=scal[i][k]; 417 imax=i; 418 kmax=k; 419 } 420 if(scal[i][k]<smin) { 421 smin=scal[i][k]; 422 imin=i; 423 kmin=k; 424 } 425 426 if(scal[i][k]>smax_value+SMALL_CONSISTENCY && grid->dzz[i][k]>DRYCELLHEIGHT) 427 maxcount++; 428 if(scal[i][k]<smin_value-SMALL_CONSISTENCY && grid->dzz[i][k]>DRYCELLHEIGHT) 429 mincount++; 430 } 431 } 432 } 433 MPI_Reduce(&mincount,&allmincount,1,MPI_INT,MPI_SUM,0,comm); 434 MPI_Reduce(&maxcount,&allmaxcount,1,MPI_INT,MPI_SUM,0,comm); 435 436 if(mincount!=0 || maxcount!=0) 437 printf("Not CWC, step: %d, proc: %d, smin = %e at i=%d,H=%e, smax = %e at i=%d,H=%e\n", 438 prop->n,myproc, 439 smin,imin,phys->h[imin]+grid->dv[imin], 440 smax,imax,phys->h[imax]+grid->dv[imax]); 441 442 if(myproc==0 && (allmincount !=0 || allmaxcount !=0)) 443 printf("Total number of CWC violations (all procs): s<s_min: %d, s>s_max: %d\n", 444 allmincount,allmaxcount); 445 } 446 }
下面介绍解读UpdateScalars函数过程：
1. 首先作为一个复杂的非静压N－S模型，变量比较多是很正常的，所以不要纠结每个变量是什么意思，能看懂就看，看不懂就猜好了。
2.要从整体入手。根据目前已知信息，这是用有限体积法求解对流扩散方程模块，而所求标量值应该就是就是第5个参数 **scal 所代表的变量。那么从函数最后一次更新 scal 变量的地方，或许能获得些许线索。
第311行：
if(grid->Nk[i]-ktop>1) TriSolve(a,b,c,d,&(scal[i][ktop]),grid->Nk[i]-ktop);
检查 TriSolve 函数的定义，原来是求解三对角方程组的解法，a，b，c 分别是系数矩阵三个对角向量，d是等号右端常向量，未知数为以 scal[i][ktop] 起始的数组。
而准备a，b，c 等系数向量时，循环变量多是按照某个三棱柱各层从上到下进行循环，所以不难看出，这个方程组求解的应该就是某个三棱柱单元体内各层标量大小。
3.
原文：http://www.cnblogs.com/li12242/p/4003350.html
