---

momentum.denis.c

Go to the documentation of this file.

/*---------------------------------------------------------------------------*/
/* SCCS Information: %W%    %G% */
/*---------------------------------------------------------------------------*/
/*  
                                  
                        o..............o
                        :              :
                    a(i-1,j)           :
                    u(i,j)             :
                 *------>------*   S11, S22(i,j), pressure, divergence
                 |      :      |   u**2, v**2, cc(i,j)
                 |      :      |       :
                 |   S12       | v(i,j):
                 ^      o..... ^.......o 
                 |             | c(i,j-1)
                 |             |
                 |             |
                 *------>------*
    
   
       SKETCH OF TOP-BOTTOM BOUNDARY CONDITIONS
    
        *-----*------*------*------*  j =  ny   p(i,ny)
        |     |      |      |      |
    ....^.....^......^......^......^....        v(i,ny) = 0
        |     |      |      |      |
        |---------------------->---
        |     |      |      |      |  j = ny-1
        |     |      |      |      ^
        |--------------------------|
        |     |      |      |      |
    
    
        |     |      |      |      |
        |--------------------------
        |     |      |      |      |
        |     |      |      |  u(n-1,2)
        |---------------------->-- |  p(n-1,2)            | Boundary conditions
        |     |      |      |      |                      |    u(i,1) = - u(i,2)
    ....^..........................^...  y = 0 v(i,2)     |    p(i,1) = p(i,2) 
        |     |      |      |      |                      |    v(i,2) = 0
        *-->--*------*------*-->-- *   j=1                   
        |     |      |      |      |   p(n-1,1)          
        ^     ^      ^      ^      ^   v(n-1,1)
        v(2,1)
    
    Sites on or outside the dotted lines are used to force boundary conditions
    There is also one dummy column on each side (not represented).


    u,v(i,j)  the physical velocities *tau/h

*/
#include "utilf.h"
#include "markers.h"
#include "extra.h"
#include "momentum.h"
#ifdef MG
  #include "mg.h"
#else
  #include "vmg.h"
#endif
#include "make_bc.h"

static void coli(real2D cc,
                 real2D nu1, real2D nu2, real2D nu3, real2D nu4, 
                 real2D c0, real2D c1, real2D c2, real2D c3, real2D c4,
                 int nx, int ny);

void momentum(real2D u, real2D v,
              real2D a, real2D c,
              real2D w1,
              real2D S11, real2D S22, real2D S12,
              real2D w2,
              real visliq,
              real gx, real gy,
              real tau,
              int nx, int ny)
{
  int i, j;
  real h;
  real gx2, gy2, tau2h;
  
  h = 1. / (nx - 2);
  tau2h = sq(tau) / h;
  gx2 = gx * tau2h;
  gy2 = gy * tau2h;
  visliq *= tau/(h*h);

  /* -----------------------------------------------------------------
     Advection term
     ------------------------------------------------------------------ */
  for (i = 2; i <= nx - 1; i++)
    {
      for (j = 2; j <= ny - 1; j++)
        if (IN(a[i-1][j]))
          w1[i][j] = 0.25*(
                      (u[i][j] + u[i][j-1])*(v[i][j] + v[i-1][j])
                    - (u[i][j+1] + u[i][j])*(v[i][j+1] + v[i-1][j+1])
                    + sq(u[i-1][j] + u[i][j]) - sq(u[i][j] + u[i+1][j])
                      /* axisymetric term */
                    - (v[i][j] + v[i][j+1] + v[i-1][j+1] + v[i-1][j])
                      *u[i][j]/((real)j - 1.5));
                     
      for (j = 3; j <= ny - 1; j++)
        if (IN(c[i][j-1]))
          w2[i][j] = 0.25*(
                      (u[i][j] + u[i][j-1])*(v[i][j] + v[i-1][j])
                    - (u[i+1][j] + u[i+1][j-1])*(v[i+1][j] + v[i][j])
                    + sq(v[i][j-1] + v[i][j]) - sq(v[i][j] + v[i][j+1])
                      /* axisymetric term */
                    - sq(2.0*v[i][j])/((real)j - 2.0));
    }

  /* ------------------------------------------------------------------
     div(S)
     Gravity term               
     ------------------------------------------------------------------ */
  for (i = 2; i <= nx - 1; i++)
    {
      for (j = 2; j <= ny - 1; j++)
        if (IN(a[i-1][j]))
          u[i][j] += w1[i][j]
            + S12[i][j+1] - S12[i][j] + S11[i][j] - S11[i-1][j]
            /* axisymetric term */
            /*
              + 0.5*(S12[i][j] + S12[i][j+1])/((real)j - 1.5)
            */
            + gx2;
      for (j = 3; j <= ny - 1; j++)
        if (IN(c[i][j-1]))
          v[i][j] += w2[i][j]
            + S12[i+1][j] - S12[i][j] + S22[i][j] - S22[i][j-1]
            /* axisymetric term */
            + (0.5*(S22[i][j] + S22[i][j-1]) - 2.*visliq*v[i][j]/
               ((real)j - 2.0))
            /((real)j - 2.0)
            + gy2;
    }
}


void pressure(real2D u, real2D v, real2D p, 
              real2D cc, real2D a, real2D c,
              real2D div,
              real2D nu1, real2D nu2, real2D nu3, real2D nu4,
              real2D c0, real2D c1, real2D c2, real2D c3, real2D c4,
#ifdef MG
              int ncycle, int npre, int npost, real pa,
#else
              real epsilon,
#endif
              int nx, int ny, int ng)
{
  int i, j;
  real h = 1.0/(nx - 2), hi2;
  real p0, p3, p4, n3, n4;
  hi2 = 1./(h*h);       

  /* ------------------------------------------------------------------
     Resolution of the pressure equation
     div=divergence
     p=pressure
     ------------------------------------------------------------------ */
  for (i = 2; i <= nx - 1; i++)
    for (j = 2; j <= ny - 1; j++)
      if (IN(cc[i][j]))
        div[i][j] = (((real)j - 1.)*v[i][j+1] - ((real)j - 2.)*v[i][j] +
          ((real)j - 1.5)*(u[i+1][j] - u[i][j]))
#ifndef MG
          *hi2
#endif
          ;
      else
        div[i][j] = 0.0;

  datbarray(div, nx, ny, "divbefore");
  datbarray(p, nx, ny, "pbefore");
  datbarray(cc, nx, ny, "cc");
  datbarray(u, nx, ny, "ubefore");
  datbarray(v, nx, ny, "vbefore");

  coli(cc,
       nu1, nu2, nu3, nu4,
       c0, c1, c2, c3, c4,
       nx, ny);
#ifdef MG
  /*
    mginit(div, p, cc, ra, rc, rei, pint, work, in, nx, ny, ng);
    for (i = 0; i < ncycle; i++)
    mgvcycle(ng, nx, ny, npre, npost);
  */
  for (i = 0; i < ncycle; i++) {
    relax(p, div, cc,
          nu1, nu2, nu3, nu4,
          c0, c1, c2, c3, c4,
          nx, ny, pa);
    bc_scalar(p, nx, ny, NULGRAD2);
  }    
#else
  mglin(div, p, cc, ra, rc, rei, epsilon, nx, ny, ng);
#endif

  datbarray(p, nx, ny, "pafter");

  /* ------------------------------------------------------------------
     Pressure correction
     ------------------------------------------------------------------ */ 
  for (i = 2; i <= nx - 1; i++) {
    for (j = 2; j <= ny - 1; j++)
      if (IN(a[i-1][j])) {
        if (IN(cc[i][j])) {
          p0 = p[i][j];
          if (IN(cc[i-1][j])) { p3 = p[i-1][j]; n3 = 1.0; }
          else { p3 = pa; n3 = -nu3[i][j]; }
        }
        else if (IN(cc[i-1][j])) {
          p3 = p[i-1][j];
          p0 = pa; n3 = nu1[i-1][j];
        }
        else {
          fprintf(stderr, "pressure(): error for u(%d,%d)\n", i, j);
          exit(1);
        }
        u[i][j] -= (p0 - p3)/n3;
      }
    for (j = 3; j <= ny - 1; j++)
      if (IN(c[i][j-1])) {
        if (IN(cc[i][j])) {
          p0 = p[i][j];
          if (IN(cc[i][j-1])) { p4 = p[i][j-1]; n4 = 1.0; }
          else { p4 = pa; n4 = -nu4[i][j]; }
        }
        else if (IN(cc[i][j-1])) {
          p4 = p[i][j-1];
          p0 = pa; n4 = nu2[i][j-1];
        }
        else {
          fprintf(stderr, "pressure(): error for v(%d,%d)\n", i, j);
          exit(1);
        }
        v[i][j] -= (p0 - p4)/n4;
      }
  }
}


void coli(real2D cc,
          real2D nu1, real2D nu2, real2D nu3, real2D nu4, 
          real2D c0, real2D c1, real2D c2, real2D c3, real2D c4,
          int nx, int ny)
{
  int i, j;
  real n1, n2, n3, n4, rj, rj12, rj_12;

  for (i = 2; i < nx; i++)
    for (j = 2; j < ny; j++) 
      if (IN(cc[i][j])) {
        n1 = IN(cc[i+1][j]) ? 1.0 : nu1[i][j];
        n2 = IN(cc[i][j+1]) ? 1.0 : nu2[i][j];
        n3 = IN(cc[i-1][j]) ? 1.0 : - nu3[i][j];
        n4 = IN(cc[i][j-1]) ? 1.0 : - nu4[i][j];

        if (n1 < 0.0 || n2 < 0.0 || n3 < 0.0 || n4 < 0.0) {
          fprintf(stderr, "coli(): error at %d,%d: n1: %g n2: %g n3: %g n4: %g\n", i, j, n1, n2, n3, n4);
        }
        
        rj = (real)j - 1.5;
        rj12 = (real)j - 1.0;
        rj_12 = (real)j - 2.0;
        c0[i][j] = n1*n2*n3*n4*(n2 + n4)/(rj*n2*n4*(n2 + n4) + 
                                          n1*n3*(rj12*n4 + rj_12*n2));
        c1[i][j] = rj/(n1*(n1 + n3));
        c2[i][j] = rj12/(n2*(n2 + n4));
        c3[i][j] = rj/(n3*(n1 + n3));
        c4[i][j] = rj_12/(n4*(n2 + n4));
      }
}


void adaptative(real2D u, real2D v, real2D divafter, real2D cc,
                int *npre, int *npost,
                real mineps, real maxeps, int nx, int ny)
{
  int i, j;
  real epsilon, h = 1.0/(nx - 2), hi2;
  hi2 = 1.0/(h*h);

  /* ------------------------------------------------------------------
     Adaptation of the multigrid for the given precision
     ------------------------------------------------------------------ */
  for (i = 2; i <= nx - 1; i++)
    for (j = 2; j <= ny - 1; j++)
      if (IN(cc[i][j]))
        divafter[i][j] = (((real)j - 1.)*v[i][j+1] - ((real)j - 2.)*v[i][j] + 
                          ((real)j - 1.5)*(u[i+1][j] - u[i][j]))*hi2;
      else
        divafter[i][j] = 0.0;
  /*
    datbarray(divafter, nx, ny, "divafter");
    */
  epsilon = rdivmax(divafter, nx, ny);
  if (epsilon > mineps) {
    (*npre)++; (*npost)++;
  }
  else if (epsilon < maxeps) {
    *npre = MAX(*npre - 1, 1);
    *npost = MAX(*npost - 1, 1);
  }
}


void adaptative_tau(real2D u, real2D v, 
                    real2D divafter, real2D cc, real2D p,
                    real *tau,
                    real mineps, real maxeps, int nx, int ny)
{
  int i, j;
  real epsilon, velmax, h = 1.0/(nx - 2), hi2;
  hi2 = 1./(h*h);

  /* ------------------------------------------------------------------
     Adaptation of the time step for the given precision
     ------------------------------------------------------------------ */
  for (i = 2; i <= nx - 1; i++)
    for (j = 2; j <= ny - 1; j++)
      if (IN(cc[i][j]))
        divafter[i][j] = (((real)j - 1.)*v[i][j+1] - ((real)j - 2.)*v[i][j] +
                          ((real)j - 1.5)*(u[i+1][j] - u[i][j]))*hi2;
      else
        divafter[i][j] = 0.0;

    datbarray(divafter, nx, ny, "divafter");

  epsilon = rdivmax(divafter, nx, ny);
  velmax = findninf(u, v, nx, ny, &i, &j);
  if (epsilon > mineps || velmax > 0.1) {
    rescale(u, v, p, 0.5, nx, ny);
    *tau *= 0.5;
  }
  else if (epsilon < maxeps && velmax < 0.05) {
    rescale(u, v, p, 2.0, nx, ny);
    *tau *= 2.;
  }
}


void rescale(real2D u, real2D v, real2D p,
             real tscale,
             int nx, int ny)
{
  int i, j;
  real tscale2 = sq(tscale);
  
  for (i = 1; i <= nx; i++)
    for (j = 1; j <= ny; j++) {
      u[i][j] *= tscale;
      v[i][j] *= tscale;
      p[i][j] *= tscale2;
    }    
}

---