momentum.c Source File

00001 /*---------------------------------------------------------------------------*/
00002 /* SCCS Information: %W%    %G% */
00003 /*---------------------------------------------------------------------------*/
00004 /*  
00005                                   
00006                         o..............o
00007                         :              :
00008                     a(i-1,j)           :
00009                     u(i,j)             :
00010                  *------>------*   S11, S22(i,j), pressure, divergence
00011                  |      :      |   u**2, v**2, cc(i,j)
00012                  |      :      |       :
00013                  |   S12       | v(i,j):
00014                  ^      o..... ^.......o 
00015                  |             | c(i,j-1)
00016                  |             |
00017                  |             |
00018                  *------>------*
00019     
00020    
00021        SKETCH OF TOP-BOTTOM BOUNDARY CONDITIONS
00022     
00023         *-----*------*------*------*  j =  ny   p(i,ny)
00024         |     |      |      |      |
00025     ....^.....^......^......^......^....        v(i,ny) = 0
00026         |     |      |      |      |
00027         |---------------------->---
00028         |     |      |      |      |  j = ny-1
00029         |     |      |      |      ^
00030         |--------------------------|
00031         |     |      |      |      |
00032     
00033     
00034         |     |      |      |      |
00035         |--------------------------
00036         |     |      |      |      |
00037         |     |      |      |  u(n-1,2)
00038         |---------------------->-- |  p(n-1,2)            | Boundary conditions
00039         |     |      |      |      |                      |    u(i,1) = - u(i,2)
00040     ....^..........................^...  y = 0 v(i,2)     |    p(i,1) = p(i,2) 
00041         |     |      |      |      |                      |    v(i,2) = 0
00042         *-->--*------*------*-->-- *   j=1                   
00043         |     |      |      |      |   p(n-1,1)          
00044         ^     ^      ^      ^      ^   v(n-1,1)
00045         v(2,1)
00046     
00047     Sites on or outside the dotted lines are used to force boundary conditions
00048     There is also one dummy column on each side (not represented).
00049 
00050 
00051     u,v(i,j)  the physical velocities *tau/h
00052 
00053 */
00054 #include "utilf.h"
00055 #include "markers.h"
00056 #include "surfaces.h"
00057 #include "extra.h"
00058 #include "momentum.h"
00059 #include "mg.h"
00060 #include "make_bc.h"
00061 #include "fluxes.h"
00062 #include "inout.h"
00063 
00064 
00065 void momentum(real2D u, real2D v,
00066               real2D cu, real2D cv, real2D ap, real2D av,
00067               real2D sxu, real2D syu, real2D sxv, real2D syv,         
00068               interface in, interface_cut cutu, interface_cut cutv,
00069               real2D rz1dr, real2D rz2dr, real2D rr1dz, real2D rr2dz,
00070               real2D Sx12,
00071               real2D S11, real2D S22, real2D S12,
00072               real2D Sy12,
00073               real visliq,
00074               real gx, real gy,
00075               real tau,
00076               int nx, int ny)
00077 {
00078   int i, j;
00079   real h, f;
00080   
00081   h = 1./(nx - 2);
00082   gx *= sq(tau)/h;
00083   gy *= sq(tau)/h;
00084   visliq *= tau/(h*h);
00085 
00086   /* -----------------------------------------------------------------
00087      Advection term
00088      ------------------------------------------------------------------ */ 
00089 #ifdef UPWIND
00090   for (i = 1; i <= nx; i++)
00091     for (j = 2; j < ny; j++) {
00092       f = 0.5*(u[i][j] + u[i+1][j]);
00093       S11[i][j] -= f > 0.0 ? u[i][j]*f : u[i+1][j]*f;
00094     }
00095   for (i = 2; i < nx; i++)
00096     for (j = 2; j <= ny; j++) {
00097       f = 0.5*(v[i][j] + v[i][j+1]);
00098       S22[i][j] -= f > 0.0 ? v[i][j]*f : v[i][j+1]*f;
00099     }
00100   for (i = 2; i <= nx; i++)
00101     for (j = 2; j <= ny; j++) {
00102       f = 0.5*(v[i][j] + v[i-1][j]);
00103       Sx12[i][j] = S12[i][j] - (f > 0.0 ? u[i][j-1]*f : u[i][j]*f);
00104       f = 0.5*(u[i][j] + u[i][j-1]);
00105       Sy12[i][j] = S12[i][j] - (f > 0.0 ? v[i-1][j]*f : v[i][j]*f);
00106     }
00107 #else
00108   for (i = 1; i <= nx; i++)
00109     for (j = 2; j < ny; j++)
00110       S11[i][j] -= 0.25*sq(u[i][j] + u[i+1][j]);
00111   for (i = 2; i < nx; i++)
00112     for (j = 2; j <= ny; j++)
00113       S22[i][j] -= 0.25*sq(v[i][j] + v[i][j+1]);
00114   for (i = 2; i <= nx; i++)
00115     for (j = 2; j <= ny; j++)
00116       Sx12[i][j] = Sy12[i][j] = S12[i][j] - 
00117         0.25*(u[i][j] + u[i][j-1])*(v[i][j] + v[i-1][j]);
00118 #endif
00119   /* ------------------------------------------------------------------
00120      div(S)
00121      Gravity term               
00122      ------------------------------------------------------------------ */
00123   interface_rzdr(in, cutu, -0.5, rz1dr, nx, ny);
00124   interface_rzdr(in, cutu, 0.5, rz2dr, nx, ny);
00125   interface_rrdz(in, cutu, -0.5, rr1dz, nx, ny);
00126   interface_rrdz(in, cutu, 0.5, rr2dz, nx, ny);
00127   for (i = 2; i <= nx; i++)
00128     for (j = 2; j < ny; j++)
00129       if (INU(i,j))
00130         u[i][j] +=
00131           ((syu[i+1][j] + rz1dr[i][j])*S11[i][j] - 
00132            (syu[i][j] + rz2dr[i][j])*S11[i-1][j] +
00133            (sxu[i][j+1] - rr1dz[i][j])*Sx12[i][j+1] - 
00134            (sxu[i][j] - rr2dz[i][j])*Sx12[i][j])/cu[i][j] +
00135           gx;
00136   interface_rzdr(in, cutv, -0.5, rz1dr, nx, ny);
00137   interface_rzdr(in, cutv, 0.5, rz2dr, nx, ny);
00138   interface_rrdz(in, cutv, -0.5, rr1dz, nx, ny);
00139   interface_rrdz(in, cutv, 0.5, rr2dz, nx, ny);
00140   for (i = 2; i < nx; i++)
00141     for (j = 3; j <= ny; j++)
00142       if (INV(i,j))
00143         v[i][j] += 
00144           ((syv[i+1][j] + rz1dr[i][j])*Sy12[i+1][j] -
00145            (syv[i][j] + rz2dr[i][j])*Sy12[i][j] + 
00146            (sxv[i][j+1] - rr1dz[i][j])*S22[i][j] -
00147            (sxv[i][j] - rr2dz[i][j])*S22[i][j-1] -
00148            av[i][j]*2.*visliq*v[i][j]/((real)j - 2.))/cv[i][j] +
00149           gy;
00150 }
00151 
00152 
00153 void pressure(real2D u, real2D v, real2D p, 
00154               real2D cu, real2D cv,
00155               real2D ap, real2D av,
00156               real2D sxp, real2D syp, 
00157               real2D syu, real2D sxv,
00158               real2D rz1dr, real2D rz2dr, real2D rr1dz, real2D rr2dz,
00159               real2D div, real2D tmp,
00160               real2D c0, real2D c1, real2D c2, real2D c3, real2D c4,
00161               interface in, real lredis,
00162               int *ncycle, int npre, int npost, real maxdiv,
00163               int nx, int ny, int ng)
00164 {
00165   int i, j;
00166 
00167   /* ------------------------------------------------------------------
00168      Resolution of the pressure equation
00169      div=divergence
00170      p=pressure
00171      ------------------------------------------------------------------ */
00172   divergence(u, v, ap, sxp, syp, rz1dr, rz2dr, rr1dz, rr2dz, div, nx, ny);
00173 #if 0
00174   datbarray(div, nx, ny, "divbefore");
00175   datbarray(u, nx, ny, "ubefore");
00176   datbarray(v, nx, ny, "vbefore");
00177   datbarray(p, nx, ny, "pbefore");
00178   datbarray(rz1dr, nx, ny, "rz1dr");
00179   datbarray(rz2dr, nx, ny, "rz2dr");
00180   datbarray(rr1dz, nx, ny, "rr1dz");
00181   datbarray(rr2dz, nx, ny, "rr2dz");
00182   datbarray(sxp, nx, ny, "sxp");
00183   datbarray(syp, nx, ny, "syp");
00184   datbarray(sxv, nx, ny, "sxv");
00185   datbarray(syu, nx, ny, "syu");
00186   datbarray(ap, nx, ny, "ap");
00187   datbarray(cv, nx, ny, "cv");
00188   datbarray(av, nx, ny, "av");
00189   datbarray(cu, nx, ny, "cu");
00190   {
00191     FILE *fptr = fopen("sxp", "wt");
00192     sx_write(sxp, 0.0, 0.0, nx, ny, fptr);
00193     fclose(fptr);
00194     
00195     fptr = fopen("sxv", "wt");
00196     sx_write(sxv, 0.0, -0.5, nx, ny, fptr);
00197     fclose(fptr);
00198 
00199     fptr = fopen("toto", "wt");
00200     interface_write_markers(in, fptr);
00201     fclose(fptr);
00202   }
00203 fill0(c0, nx, ny);
00204 fill0(c1, nx, ny);
00205 fill0(c2, nx, ny);
00206 fill0(c3, nx, ny);
00207 fill0(c4, nx, ny);
00208 #endif
00209   coli(sxp, syp, syu, sxv, cu, cv, ap, av, 
00210        rz1dr, rz2dr, rr1dz, rr2dz,
00211        c0, c1, c2, c3, c4, nx, ny, 1.0);
00212 #if 0
00213   datbarray(c0, nx, ny, "c0");
00214   datbarray(c1, nx, ny, "c1");
00215   datbarray(c2, nx, ny, "c2");
00216   datbarray(c3, nx, ny, "c3");
00217   datbarray(c4, nx, ny, "c4");
00218 #endif
00219  
00220 #if 1
00221    copy(tmp, p, nx, ny);
00222   *ncycle = mgsolve(div, p, ap,
00223                     c0, c1, c2, c3, c4, 
00224                     in, lredis, maxdiv, nx, ny, ng);
00225 #else
00226   for (i = 0; i < 20; i++) {
00227     relax(p, div, ap, c0, c1, c2, c3, c4, nx, ny);
00228     bc_scalar(p, nx, ny, NULGRAD2);
00229   }
00230 #endif
00231 
00232   /* ------------------------------------------------------------------
00233      Pressure correction
00234      ------------------------------------------------------------------ */ 
00235   for (i = 2; i <= nx; i++)
00236     for (j = 2; j < ny; j++)
00237       if (INU(i,j))
00238         u[i][j] += (syu[i][j]*p[i-1][j] - syu[i+1][j]*p[i][j])/cu[i][j];
00239   for (i = 2; i < nx; i++)
00240     for (j = 3; j < ny; j++) /* !!!! <= ny when No top wall !!!! */
00241       if (INV(i,j))
00242         v[i][j] += ((sxv[i][j] + av[i][j]/2.)*p[i][j-1] - 
00243                     (sxv[i][j+1] - av[i][j]/2.)*p[i][j])/cv[i][j];
00244 #if 0
00245     datbarray(u, nx, ny, "uafter");
00246     datbarray(v, nx, ny, "vafter");
00247     datbarray(p, nx, ny, "pafter");
00248 #endif
00249 }
00250 
00251 
00252 void coli(real2D sxp, real2D syp,
00253           real2D syu, real2D sxv,
00254           real2D cu, real2D cv, real2D ap, real2D av,
00255           real2D rz1dr, real2D rz2dr,
00256           real2D rr1dz, real2D rr2dz,
00257           real2D c0, real2D c1, real2D c2, real2D c3, real2D c4,
00258           int nx, int ny, real scale)
00259 {
00260   int i, j;
00261   real syp1, syp3, sxp2, sxp4;
00262 
00263   for (i = 2; i < nx; i++)
00264     for (j = 2; j < ny; j++) 
00265       if (INP(i,j)) {
00266         syp1 = scale*(syp[i+1][j] + rz1dr[i][j])/cu[i+1][j];
00267         syp3 = scale*(syp[i][j] + rz2dr[i][j])/cu[i][j];
00268 
00269         c1[i][j] = syu[i+2][j]*syp1;
00270         c3[i][j] = syu[i][j]*syp3;
00271 
00272         sxp2 = j == ny-1 ? 0.0 : scale*(sxp[i][j+1] - rr1dz[i][j])/cv[i][j+1];
00273         sxp4 = j == 2 ? 0.0 : scale*(sxp[i][j] - rr2dz[i][j])/cv[i][j];
00274 
00275         c2[i][j] = (sxv[i][j+2] - 0.5*av[i][j+1])*sxp2;
00276         c4[i][j] = (sxv[i][j] + 0.5*av[i][j])*sxp4;
00277         
00278         c0[i][j] = syu[i+1][j]*(syp1 + syp3) + sxv[i][j+1]*(sxp2 + sxp4) +
00279           0.5*(av[i][j+1]*sxp2 - av[i][j]*sxp4);
00280 #ifdef CHECK_DOMINANCE
00281         if (fabs(c0[i][j]) < fabs(c1[i][j])*INP(i+1,j) + 
00282             fabs(c2[i][j])*INP(i,j+1) + 
00283             fabs(c3[i][j])*INP(i-1,j) + 
00284             fabs(c4[i][j])*INP(i,j-1)) {
00285           printf("WARNING: coli(%d,%d): c0: %g c1: %g c2: %g c3: %g c4: %g\n",
00286                  i, j, c0[i][j], c1[i][j]*INP(i+1,j), c2[i][j]*INP(i,j+1), 
00287                  c3[i][j]*INP(i-1,j), c4[i][j]*INP(i,j-1));
00288           c0[i][j] = fabs(c1[i][j]) + fabs(c2[i][j]) + 
00289             fabs(c3[i][j]) + fabs(c4[i][j]);
00290           printf("         setting c0 to %g\n", c0[i][j]);
00291         }
00292 #endif
00293 #if 0
00294         if (i == 150 && j == 3 && scale == 1.0) {
00295           printf("syp[%d][%d]: %g rz1dr[%d][%d]: %g cu[%d][%d]: %g\n",
00296                  i+1, j, syp[i+1][j], i, j, rz1dr[i][j], i+1, j, cu[i+1][j]);
00297           printf("syp[%d][%d]: %g rz2dr[%d][%d]: %g cu[%d][%d]: %g\n",
00298                  i, j, syp[i][j], i, j, rz2dr[i][j], i, j, cu[i][j]);
00299           printf("syu[%d][%d]: %g c1[%d][%d]: %g\n", 
00300                  i+2, j, syu[i+2][j], i, j, c1[i][j]);
00301           printf("syu[%d][%d]: %g c3[%d][%d]: %g\n", 
00302                  i, j, syu[i][j], i, j, c3[i][j]);
00303           printf("\n");
00304           printf("sxp[%d][%d]: %g rr1dz[%d][%d]: %g cv[%d][%d]: %g\n",
00305                  i, j+1, sxp[i][j+1], i, j, rr1dz[i][j], i, j+1, cv[i][j+1]);
00306           printf("sxp[%d][%d]: %g rr2dz[%d][%d]: %g cv[%d][%d]: %g\n",
00307                  i, j, sxp[i][j], i, j, rr2dz[i][j], i, j, cv[i][j]);
00308           printf("sxv[%d][%d]: %g av[%d][%d]: %g c2[%d][%d]: %g\n", 
00309                  i, j+2, sxv[i][j+2], i, j+1, av[i][j+1], i, j, c2[i][j]);
00310           printf("sxv[%d][%d]: %g av[%d][%d]: %g c4[%d][%d]: %g\n", 
00311                  i, j, sxv[i][j], i, j, av[i][j], i, j, c4[i][j]);
00312           printf("\n");
00313           printf("syu[%d][%d]: %g sxv[%d][%d]: %g c0[%d][%d]: %g\n",
00314                  i+1, j, syu[i+1][j], i, j+1, sxv[i][j+1], i, j, c0[i][j]);
00315         } 
00316 #endif
00317       }
00318 }
00319 
00320 
00321 void adaptative(real2D u, real2D v, real2D p,
00322                 real *tau,
00323                 int ncycle, int nx, int ny)
00324 {
00325   real umax, vmax;
00326   int i, j;
00327 
00328   umax = fmodmax(u, nx, ny, &i, &j);
00329   vmax = fmodmax(v, nx, ny, &i, &j);
00330   if (vmax > umax) umax = vmax;
00331   if (ncycle > MAXCYCLE || umax > MAXVEL) {
00332     rescale(u, v, p, 1./DIVIDE, nx, ny);
00333     *tau /= DIVIDE;
00334   }
00335   else if (ncycle < 2 && umax*MULT < MAXVEL) {
00336     rescale(u, v, p, MULT, nx, ny);
00337     *tau *= MULT;
00338   }
00339 }
00340 
00341 
00342 void rescale(real2D u, real2D v, real2D p,
00343              real tscale,
00344              int nx, int ny)
00345 {
00346   int i, j;
00347   real tscale2 = sq(tscale);
00348   
00349   for (i = 1; i <= nx; i++)
00350     for (j = 1; j <= ny; j++) {
00351       if (u[i][j] != UNDEFINED) u[i][j] *= tscale;
00352       if (v[i][j] != UNDEFINED) v[i][j] *= tscale;
00353       p[i][j] *= tscale2;
00354     }
00355 }
00356 
00357 
00358 void divergence(real2D u, real2D v, real2D ap,
00359                 real2D sxp, real2D syp, 
00360                 real2D rz1dr, real2D rz2dr,
00361                 real2D rr1dz, real2D rr2dz,
00362                 real2D div, int nx, int ny)
00363 {
00364   int i, j;
00365 
00366   for (i = 2; i < nx; i++)
00367     for (j = 2; j < ny; j++)
00368       if (INP(i,j)) {
00369         div[i][j] = (syp[i+1][j] + rz1dr[i][j])*u[i+1][j] - 
00370           (syp[i][j] + rz2dr[i][j])*u[i][j] + 
00371           (sxp[i][j+1] - rr1dz[i][j])*v[i][j+1] - 
00372           (sxp[i][j] - rr2dz[i][j])*v[i][j];
00373       }
00374       else
00375         div[i][j] = 0.0;
00376 }
00377 
00378 
00379 
00380 
00381 
00382