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 tridiagonal117     // 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 schemes124           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 advection134       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 advection139       b[(grid->Nk[i]-1)-ktop]+=c[(grid->Nk[i]-1)-ktop];140 141     // Implicit vertical diffusion terms142       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 layer156       if(ktop<grid->Nk[i]-1) {157       // Top cell diffusion158           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 diffusion162           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 tridiagonal187     // that use the new velocity188       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 schemes194           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 diffusion198       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 cell208           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 cell216           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       else244           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 step248       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               else264                   sp=phys->stmp[nc1];265           } else266               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                   else278                       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 < -d297     /*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           else317               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.h328   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   }