SUNTANS 及 FVCOM 对流扩散方程求解简介［TBC］ |
您所在的位置：网站首页 › 对流扩散系数单位是什么 › SUNTANS 及 FVCOM 对流扩散方程求解简介［TBC］ |
SUNTANS 及 FVCOM 对流扩散方程求解简介［TBC］
|
最近接到一个任务，就是解决FVCOM中对流扩散计算不守衡问题。导师认为是其求解时候水平和垂向计算分开求解所导致的，目前我也没搞清到底有什么问题，反正就是让把SUNTANS的对流扩散计算挪到FVCOM中，下面就把 SUNTANS 和 FVCOM 数值求解的过程贴出来，备忘 SUNTANS模型求解过程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;iNc;i++)
64 for(k=0;kNk[i];k++)
65 Cn[i][k]=0;
66 } else
67 fab=1.5;
68
69 for(i=0;istmp[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];iptrcelldist[1];iptr++) {
76 for(iptr=grid->celldist[0];iptrcelldist[2];iptr++) {
77 i = grid->cellp[iptr];
78
79 for(k=grid->ctop[i];kNk[i];k++)
80 phys->stmp2[i][k]=0;
81 }
82
83 if(boundary_scal) {
84 for(jptr=grid->edgedist[2];jptredgedist[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];kNk[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];iptrcelldist[1];iptr++) {
102 for(iptr=grid->celldist[0];iptrcelldist[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];kctopold[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;kNk[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;kNk[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;kNk[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;kNk[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;kNk[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(ktopNk[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;kNk[i];k++)
168 d[k-ktop]=grid->dzzold[i][k]*phys->stmp[i][k];
169 if(src1)
170 for(k=ktop+1;kNk[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]ctop[i]) {
175 for(k=grid->ctopold[i];kctop[i];k++)
176 d[0]+=grid->dzzold[i][k]*phys->stmp[i][k];
177 if(src1)
178 for(k=grid->ctopold[i];kctop[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;kNk[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;kNk[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(ktopNk[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]ctopold[i]) {
226 for(k=grid->ctop[i];kctopold[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;kNk[i];k++)
229 d[k-grid->ctopold[i]]+=(1-fab)*Cn[i][k];
230 } else {
231 for(k=grid->ctopold[i];kctop[i];k++)
232 d[0]+=(1-fab)*Cn[i][k];
233 for(k=grid->ctop[i]+1;kNk[i];k++)
234 d[k-grid->ctop[i]]+=(1-fab)*Cn[i][k];
235 }
236
237 for(k=0;kctop[i];k++)
238 Cn[i][k]=0;
239
240 if(src2)
241 for(k=grid->ctop[i];kNk[i];k++)
242 Cn[i][k-ktop]=dt*src2[i][k]*grid->dzzold[i][k];
243 else
244 for(k=grid->ctop[i];kNk[i];k++)
245 Cn[i][k]=0;
246
247 // Now create the source term for the current time step
248 for(k=0;kNk[i];k++)
249 ap[k]=0;
250
251 for(nf=0;nfnfaces[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;kNke[ne];k++)
270 temp[k]=UpWind(phys->utmp2[ne][k],
271 phys->stmp[nc1][k],
272 sp[k]);
273 } else {
274 for(k=0;kNke[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;kNke[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;kNk[i];k++)
287 Cn[i][k-ktop]-=ap[k];
288
289 for(k=0;kctop[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;kNk[i];k++)
303 ap[k]=Cn[i][k-ktop];
304 for(k=0;kctop[i];kctop[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;kctop[i];k++)
321 scal[i][k]=0;
322
323 for(k=grid->ctop[i];kctopold[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;iNc;i++) {
334 for(k=grid->ctop[i];kNk[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]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];iptrcelldist[1];iptr++) {
357 for(iptr=grid->celldist[0];iptrcelldist[2];iptr++) {
358 i = grid->cellp[iptr];
359
360 flag=0;
361 for(nf=0;nfnfaces[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;kNk[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;nfnfaces[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];iptrcelldist[1];iptr++) {
402 for(iptr=grid->celldist[0];iptrcelldist[2];iptr++) {
403 i = grid->cellp[iptr];
404
405 flag=0;
406 for(nf=0;nfnfaces[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];kNk[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]smax_value+SMALL_CONSISTENCY && grid->dzz[i][k]>DRYCELLHEIGHT)
427 maxcount++;
428 if(scal[i][k]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): ss_max: %d\n",
444 allmincount,allmaxcount);
445 }
446 }
View Code
下面介绍解读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. 将程序数值离散过程和理论结合起来，了解程序细节 FVCOM 模型求解过程
FVCOM 也是使用有限体积方法，但是求解和 SUNTANS 有很大不同。由于FVCOM并没有介绍对流扩散方程求解具体过程的文献，这时程序看起来就比较头疼，只能全部通过读程序来一点点理解。 FVCOM 扩散方程计算主要在程序 mod_scal.F 中，模块内全部程序如下
1 !/===========================================================================/
2 ! Copyright (c) 2007, The University of Massachusetts Dartmouth
3 ! Produced at the School of Marine Science & Technology
4 ! Marine Ecosystem Dynamics Modeling group
5 ! All rights reserved.
6 !
7 ! FVCOM has been developed by the joint UMASSD-WHOI research team. For
8 ! details of authorship and attribution of credit please see the FVCOM
9 ! technical manual or contact the MEDM group.
10 !
11 !
12 ! This file is part of FVCOM. For details, see http://fvcom.smast.umassd.edu
13 ! The full copyright notice is contained in the file COPYRIGHT located in the
14 ! root directory of the FVCOM code. This original header must be maintained
15 ! in all distributed versions.
16 !
17 ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
18 ! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
19 ! THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
20 ! PURPOSE ARE DISCLAIMED.
21 !
22 !/---------------------------------------------------------------------------/
23 ! CVS VERSION INFORMATION
24 ! $Id$
25 ! $Name$
26 ! $Revision$
27 !/===========================================================================/
28
29 !=======================================================================
30 ! FVCOM Scalar Module
31 !
32 ! contains methods:
33 ! Adv_Scal => Advect a Scalar Quantity
34 ! Vdif_Scal => Vertical Diffusion of Scalar Quantity
35 ! Bcond_Scal_OBC => Open Boundary Condition for Scalar
36 ! Bcond_Scal_PTsource => Point Sources of Scalar
37 !=======================================================================
38 Module Scalar
39
40 logical, parameter :: debug = .true.
41
42 contains
43 !==============================================================================|
44 ! Calculate Horizontal Advection and Diffusion For Scalar (f) |
45 !==============================================================================|
46 Subroutine Adv_Scal(f,fn,d_fdis,fdis,d_fflux,fflux_obc,deltat,source)
47 !------------------------------------------------------------------------------|
48
49 use all_vars
50 use lims, only: m,mt,n,nt,kbm1,kb
51 use bcs
52 use mod_obcs
53 # if defined (MULTIPROCESSOR)
54 use mod_par
55 # endif
56 # if defined (WET_DRY)
57 use mod_wd
58 # endif
59
60 # if defined (THIN_DAM)
61 use mod_dam, only : kdam,N_DAM_MATCH,IS_DAM
62 # endif
63
64 implicit none
65 real(sp), intent(in ), dimension(0:mt,kb) :: f
66 real(sp), intent(out), dimension(0:mt,kb) :: fn
67 integer , intent(in ) :: d_fdis
68 real(sp), intent(in ), dimension(d_fdis) :: fdis
69 integer , intent(in ) :: d_fflux
70 real(sp), intent(out), dimension(d_fflux,kbm1) :: fflux_obc
71 real(sp), intent(in ) :: deltat
72 logical , intent(in ) :: source
73
74 !----------------local--------------------------------------
75 real(sp), dimension(0:mt,kb) :: xflux,xflux_adv
76 real(sp), dimension(m) :: pupx,pupy,pvpx,pvpy
77 real(sp), dimension(m) :: pfpx,pfpy,pfpxd,pfpyd,viscoff
78 real(sp), dimension(3*nt) :: dtij
79 real(sp), dimension(3*nt,kbm1) :: uvn
80 real(sp), dimension(kb) :: vflux
81 real(sp) :: utmp,vtmp,sitai,ffd,ff1,x11,y11,x22,y22,x33,y33
82 real(sp) :: tmp1,tmp2,xi,yi
83 real(sp) :: dxa,dya,dxb,dyb,fij1,fij2,un
84 real(sp) :: txx,tyy,fxx,fyy,viscof,exflux,temp,fpoint
85 real(sp) :: fact,fm1,fmean
86 integer :: i,i1,i2,ia,ib,j,j1,j2,k,jtmp,jj
87 # if defined (SPHERICAL)
88 real(sp) :: ty,txpi,typi
89 # endif
90
91 # if defined (THIN_DAM)
92 INTEGER :: NX
93 real(sp) :: tmpflx
94 real(sp),dimension(kb) :: wvel
95 # endif
96
97
98 !------------------------------------------------------------------------------!
99
100 !-------------------------------------------------------
101 !Calculate Mean Values
102 !-------------------------------------------------------
103
104 fmean = sum(f(1:m,1:kbm1))/float(m*kbm1)
105
106 !-------------------------------------------------------
107 !Initialize Multipliers to Control Horizontal Diff
108 !-------------------------------------------------------
109
110 fact = 0.0_sp
111 fm1 = 1.0_sp
112 if(HORIZONTAL_MIXING_TYPE == 'closure') then
113 fact = 1.0_sp
114 fm1 = 0.0_sp
115 end if
116
117 !-------------------------------------------------------
118 !Initialize Fluxes
119 !-------------------------------------------------------
120 xflux = 0.0_sp
121 xflux_adv = 0.0_sp
122
123 !-------------------------------------------------------
124 !Calculate Normal Velocity on Control Volume Edges
125 !-------------------------------------------------------
126 !!# if !defined (WET_DRY)
127 do i=1,ncv
128 i1=ntrg(i)
129 dtij(i)=dt1(i1)
130 do k=1,kbm1
131 uvn(i,k) = v(i1,k)*dltxe(i) - u(i1,k)*dltye(i)
132
133 # if defined(PLBC)
134 uvn(i,k) = - u(i1,k)*dltye(i)
135 # endif
136
137 end do
138 end do
139 !!# else
140 !! do i=1,ncv
141 !! i1=ntrg(i)
142 !! dtij(i)=dt1(i1)
143 !! do k=1,kbm1
144 !! uvn(i,k) = vs(i1,k)*dltxe(i) - us(i1,k)*dltye(i)
145 !! end do
146 !! end do
147 !!# endif
148
149 !
150 !--Calculate the Advection and Horizontal Diffusion Terms----------------------!
151 !
152
153 do k=1,kbm1
154 pfpx = 0.0_sp
155 pfpy = 0.0_sp
156 pfpxd = 0.0_sp
157 pfpyd = 0.0_sp
158 do i=1,m
159 do j=1,ntsn(i)-1
160 i1=nbsn(i,j)
161 i2=nbsn(i,j+1)
162
163 # if defined (WET_DRY)
164 IF(ISWETN(I1) == 0 .AND. ISWETN(I2) == 1)THEN
165 FFD=0.5_SP*(f(I,K)+f(I2,K))
166 FF1=0.5_SP*(f(I,K)+f(I2,K))
167 ELSE IF(ISWETN(I1) == 1 .AND. ISWETN(I2) == 0)THEN
168 FFD=0.5_SP*(f(I1,K)+f(I,K))
169 FF1=0.5_SP*(f(I1,K)+f(I,K))
170 ELSE IF(ISWETN(I1) == 0 .AND. ISWETN(I2) == 0)THEN
171 FFD=0.5_SP*(f(I,K)+f(I,K))
172 FF1=0.5_SP*(f(I,K)+f(I,K))
173 ELSE
174 FFD=0.5_SP*(f(I1,K)+f(I2,K))
175 FF1=0.5_SP*(f(I1,K)+f(I2,K))
176 END IF
177 # else
178 ffd=0.5_sp*(f(i1,k)+f(i2,k)) !-fmean1(i1,k)-fmean1(i2,k))
179 ff1=0.5_sp*(f(i1,k)+f(i2,k))
180 # endif
181
182 # if defined (SPHERICAL)
183 ty=0.5_sp*(vy(i1)+vy(i2))
184 txpi=(vx(i2)-vx(i1))*tpi*cos(deg2rad*ty)
185 typi=(vy(i1)-vy(i2))*tpi
186 pfpx(i)=pfpx(i)+ff1*typi
187 pfpy(i)=pfpy(i)+ff1*txpi
188 pfpxd(i)=pfpxd(i)+ffd*typi
189 pfpyd(i)=pfpyd(i)+ffd*txpi
190 # else
191 pfpx(i) = pfpx(i) +ff1*(vy(i1)-vy(i2))
192 pfpy(i) = pfpy(i) +ff1*(vx(i2)-vx(i1))
193 pfpxd(i)= pfpxd(i)+ffd*(vy(i1)-vy(i2))
194 pfpyd(i)= pfpyd(i)+ffd*(vx(i2)-vx(i1))
195 # endif
196 end do
197
198 ! gather all neighboring control volumes connecting at dam node
199 # if defined (THIN_DAM)
200 IF(IS_DAM(I)==1.AND.K 0) then
812 do i=1,iobcn
813 j=i_obc_n(i)
814 j1=next_obc(i)
815 f2d=0.0_sp
816 f2d_next=0.0_sp
817 xflux2d=0.0_sp
818 do k=1,kbm1
819 f2d=f2d+f(j,k)*dz(j,k)
820 f2d_next=f2d_next+fn(j1,k)*dz(j1,k)
821 xflux2d=xflux2d+fflux_obc(i,k)*dz(j,k)
822 end do
823
824 if(uard_obcn(i) > 0.0_sp) then
825 tmp=xflux2d+f2d*uard_obcn(i)
826 f2d_obc=(f2d*dt(j)-tmp*deltat/art1(j))/d(j)
827 do k=1,kbm1
828 fn(j,k)=fn(j1,k) !f2d_obc+(fn(j1,k)-f2d_next)
829 end do
830 else
831 do k=1,kbm1
832 fn(j,k) = f(j,k)-alpha_nudge*(f(j,k)-f_obc(i))
833 end do
834 end if
835 end do
836 endif
837
838 return
839 End Subroutine Bcond_Scal_OBC
840 !==============================================================================|
841 !==============================================================================|
842
843 Subroutine fct_sed(f,fn)
844 !==============================================================================|
845 USE ALL_VARS
846 USE MOD_UTILS
847 USE BCS
848 USE MOD_OBCS
849 IMPLICIT NONE
850 real(sp), intent(inout), dimension(0:mt,kb) :: fn
851 real(sp), intent(in), dimension(0:mt,kb) :: f
852 REAL(SP):: SMAX,SMIN
853 INTEGER :: I,J,K,K1
854 !==============================================================================|
855 IF(DBG_SET(DBG_SBR)) WRITE(IPT,*)"Start: fct_sed"
856
857 nodes: DO I=1,M
858
859 ! SKIP OPEN BOUNDARY NODES
860 IF(IOBCN > 0)THEN
861 DO J=1,IOBCN
862 IF(I == I_OBC_N(J)) CYCLE nodes
863 END DO
864 END IF
865
866 ! SKIP RIVER INFLOW POINTS
867 IF(NUMQBC > 0)THEN
868 DO J=1,NUMQBC
869 IF(RIVER_INFLOW_LOCATION == 'node')THEN
870 IF(I == INODEQ(J)) CYCLE nodes
871 END IF
872 IF(RIVER_INFLOW_LOCATION == 'edge')THEN
873 IF(I == N_ICELLQ(J,1) .OR. I == N_ICELLQ(J,2)) CYCLE nodes
874 END IF
875 END DO
876 END IF
877
878 ! SKIP GROUND WATER INFLOW POINTS
879 IF(BFWDIS(I) .GT. 0.0_SP .and. GROUNDWATER_SALT_ON) CYCLE nodes
880
881 K1 = 1
882 IF(PRECIPITATION_ON) K1 = 2
883 ! DO K=1,KBM1
884 DO K=K1,KBM1
885 SMAX = MAXVAL(f(NBSN(I,1:NTSN(I)),K))
886 SMIN = MINVAL(f(NBSN(I,1:NTSN(I)),K))
887
888 IF(K == 1)THEN
889 SMAX = MAX(SMAX,(f(I,K)*DZ(I,K+1)+f(I,K+1)*DZ(I,K))/ &
890 (DZ(I,K)+DZ(I,K+1)))
891 SMIN = MIN(SMIN,(f(I,K)*DZ(I,K+1)+f(I,K+1)*DZ(I,K))/ &
892 (DZ(I,K)+DZ(I,K+1)))
893 ELSE IF(K == KBM1)THEN
894 SMAX = MAX(SMAX,(f(I,K)*DZ(I,K-1)+f(I,K-1)*DZ(I,K))/ &
895 (DZ(I,K)+DZ(I,K-1)))
896 SMIN = MIN(SMIN,(f(I,K)*DZ(I,K-1)+f(I,K-1)*DZ(I,K))/ &
897 (DZ(I,K)+DZ(I,K-1)))
898 ELSE
899 SMAX = MAX(SMAX,(f(I,K)*DZ(I,K-1)+f(I,K-1)*DZ(I,K))/ &
900 (DZ(I,K)+DZ(I,K-1)), &
901 (f(I,K)*DZ(I,K+1)+f(I,K+1)*DZ(I,K))/ &
902 (DZ(I,K)+DZ(I,K+1)))
903 SMIN = MIN(SMIN,(f(I,K)*DZ(I,K-1)+f(I,K-1)*DZ(I,K))/ &
904 (DZ(I,K)+DZ(I,K-1)), &
905 (f(I,K)*DZ(I,K+1)+f(I,K+1)*DZ(I,K))/ &
906 (DZ(I,K)+DZ(I,K+1)))
907 END IF
908
909 IF(SMIN-fn(I,K) > 0.0_SP)fn(I,K) = SMIN
910 IF(fn(I,K)-SMAX > 0.0_SP)fn(I,K) = SMAX
911
912 END DO
913 END DO nodes
914
915 WHERE(fn < 0.0_SP)fn=0.0_SP
916
917 IF(DBG_SET(DBG_SBR)) WRITE(IPT,*)"End: fct_sed"
918 End Subroutine fct_sed
919
920 !==========================================================================
921 ! Calculate Fluxes for Vertical Advection Equation
922 ! n: number of cells
923 ! c: scalar variable (1:n)
924 ! w: velocity field at cell interfaces (1:n+1)
925 ! note: we dont use face normals to construct inflow/outflow
926 ! thus we add dissipation term instead of subtracting because
927 ! positive velocity is up while computational coordinates increase
928 ! down towards bottom.
929 !==========================================================================
930 Subroutine Calc_VFlux(n,c,w,flux)
931 use mod_prec
932 implicit none
933 integer , intent(in ) :: n
934 real(sp), intent(in ) :: c(n)
935 real(sp), intent(in ) :: w(n+1)
936 real(sp), intent(out) :: flux(n+1)
937 real(sp) :: conv(n+1),diss(n+1)
938 real(sp) :: cin(-1:n+2)
939 real(sp) :: dis4
940 integer :: i
941
942 !transfer to working array
943 cin(1:n) = c(1:n)
944
945 !surface bcs (no flux)
946 cin(0) = -cin(1)
947 cin(-1) = -cin(2)
948
949 !bottom bcs (no flux)
950 cin(n+1) = -cin(n)
951 cin(n+2) = -cin(n-1)
952
953 !flux computation
954 do i=1,n+1
955 dis4 = .5*abs(w(i))
956 conv(i) = w(i)*(cin(i)+cin(i-1))/2.
957 diss(i) = dis4*(cin(i)-cin(i-1)-lim(cin(i+1)-cin(i),cin(i-1)-cin(i-2)))
958 flux(i) = conv(i)+diss(i)
959 end do
960
961 End Subroutine Calc_VFlux
962
963 !==========================================================================
964 ! Calculate LED Limiter L(u,v)
965 !==========================================================================
966 Function Lim(a,b)
967 use mod_prec
968 real(sp) lim,a,b
969 real(sp) q,R
970 real(sp) eps
971 eps = epsilon(eps)
972
973 q = 0. !1st order
974 q = 1. !minmod
975 q = 2. !van leer
976
977 R = abs( (a-b)/(abs(a)+abs(b)+eps) )**q
978 lim = .5*(1-R)*(a+b)
979
980 End Function Lim
981
982
983 End Module Scalar
View Code
为了更快速的了解计算过程，自己设置了一个只有7个节点、6个单元的简单地形，如下图所示，然后通过 printf 的方法快速了解每个变量含义。
首先介绍 FVCOM 控制方程离散，虽然其也是按照有限体积方法思想将积分降维，但是只是在平面二维方向上使用有限体积法，垂向计算使用的是有限差分法。 控制方程
控制方程离散过程
最终更新标量值Ci的程序为： fn(i,k)=(f(i,k)-xflux(i,k)/art1(i)*(deltat/dt(i)))*(dt(i)/dtfa(i)) 这里deltat表示时间步长，dt为上个时间步水深，dtfa为新计算水深，这里之所以需要除以水深是因为垂向梯度项计算需要。 数值求解（结合程序分析）首先说一下 FVCOM 选取控制体的方法，如下图所示。例如节点1所在控制体，就是由下边4条红线和上边2条黑线所组成的6面体，而节点6则是由周围12条红线组成。
FVCOM计算时候也不是按照控制体进行循环，而是按照控制边的个数循环，也就是说，像节点1和节点6这种相邻控制体，两条红色邻边各只计算一次，控制边的通量在节点1和6控制体内是大小相同、符号相反的。 1. 水平对流项计算
这里前面一项 -un*dtij(i)*((1.0_sp+sign(1.0_sp,un))*fij2+(1.0_sp-sign(1.0_sp,un))*fij1)*0.5_sp 毫无疑问就是水平对流项，((1.0_sp+sign(1.0_sp,un))*fij2+(1.0_sp-sign(1.0_sp,un))*fij1)*0.5_sp 表示其采用的是迎风格式，控制边上 fij1 及 fij2 计算采用泰勒公式达到2阶精度，泰勒公式中的一阶偏导计算在后面统一介绍
txx=0.5_sp*(pfpxd(ia)+pfpxd(ib))*viscof tyy=0.5_sp*(pfpyd(ia)+pfpyd(ib))*viscof fxx=-dtij(i)*txx*dltye(i) fyy= dtij(i)*tyy*dltxe(i) exflux=-un*dtij(i)* & ((1.0_sp+sign(1.0_sp,un))*fij2+(1.0_sp-sign(1.0_sp,un))*fij1)*0.5_sp+fxx+fyy 其中 viscof 为水平扩散系数，一阶偏导采用控制边相邻两个节点平均值。 fxx 和 fyy 中还包含了水深dtij，后面会在计算流量时除去。 exflux=-un*dtij(i)* & ((1.0_sp+sign(1.0_sp,un))*fij2+(1.0_sp-sign(1.0_sp,un))*fij1)*0.5_sp+fxx+fyy xflux(ia,k)=xflux(ia,k)+exflux xflux(ib,k)=xflux(ib,k)-exflux fn(i,k)=(f(i,k)-xflux(i,k)/art1(i)*(deltat/dt(i)))*(dt(i)/dtfa(i)) 3. 一阶偏导项计算一阶偏导计算也采用格林公式将面积分降维化为线积分计算，但是平面积分所采用的控制体和上面不同，这里采用和节点相邻的所有三角形单元计算，对于边界点来说，只取相邻的两个单元。
对应的程序如下 pfpx(i) = pfpx(i) +ff1*(vy(i1)-vy(i2)) pfpy(i) = pfpy(i) +ff1*(vx(i2)-vx(i1)) pfpxd(i)= pfpxd(i)+ffd*(vy(i1)-vy(i2)) pfpyd(i)= pfpyd(i)+ffd*(vx(i2)-vx(i1)) …… PFPX(I)=PFPX(I)/ART2(I) PFPY(I)=PFPY(I)/ART2(I) PFPXD(I)=PFPXD(I)/ART2(I) PFPYD(I)=PFPYD(I)/ART2(I)对于水平对流项中控制边上Ci的二阶精度计算，只需按照泰勒公式计算即可 fij1=f(ia,k)+dxa*pfpx(ia)+dya*pfpy(ia) fij2=f(ib,k)+dxb*pfpx(ib)+dyb*pfpy(ib)
|
【本文地址】
公司简介
联系我们
今日新闻 |
点击排行 |
|
推荐新闻 |
图片新闻 |
|
专题文章 |