---

make_bc_wall.c

Go to the documentation of this file.

/*---------------------------------------------------------------------------*/
/* SCCS Information: @(#)make_bc.c      1.1    3/7/96 */
/*---------------------------------------------------------------------------*/
/* Implements the following boundary conditions
  
   periodic in x direction. 
   slip velocity           on y=0,1
   V.C = 0                 on y=0,1 
 
*/   
#include "utilf.h"
#include "stypes.h"
#include "markers.h"
#include "extra.h"
#include "extra1.h"
#include "make_bc.h"
#include "make_bc_periodic.h"
#include "inout.h"

#define WALL_LEFT
/*  #define WALL_RIGHT */
#define WALL_TOP


static real rp_pressure(real pe, real R, real x, real y, int nx, int ny)
{
  real r = sqrt(sq(x - 0.5*nx - 1.0) + sq(y - 2.0));
#if 1
  return pe;
#else
  return pe*(1. - R/r);
#endif
}


void bc_vector_bound(real2D u, real2D v, int nx, int ny)
{
  int i, j;
  /* top and bottom walls */
  for (i = 1; i <= nx; i++) {
    u[i][1] = u[i][2]; /* axisymmetric */
    v[i][2] = 0.0;     /* axisymmetric */
    u[i][ny] = u[i][ny-1];
#ifdef WALL_TOP
    v[i][ny] = 0.0;
#else
    v[i][ny+1] = v[i][ny];
#endif
  }
  /* left and right walls */
  for (j = 1; j <= ny; j++) {
#ifdef WALL_LEFT
    u[2][j] = 0.0;
#else
    u[1][j] = u[2][j];
#endif
#ifdef WALL_RIGHT
    u[nx][j] = 0.0;
#else
    u[nx+1][j] = u[nx][j];
#endif
    v[1][j] = v[2][j];
    v[nx][j] = v[nx-1][j];
  }
}


void bc_vector_div(real2D u, real2D v, real2D ap,
                   interface in,
                   int nx, int ny)
{
  int i, j;
  real eu, ev;

  /* extrapolate the interfacial cells */
  for (i = 3; i < nx; i++) {
    for (j = 2; j < ny; j++)
      if (!INU(i,j) && BU(i,j)) {
        /*
        if (INU(i+1,j) && INU(i-1,j)) {
          u[i][j] = (u[i+1][j] + u[i-1][j])/2.;
          printf("average u: %d,%d + %d,%d\n", i+1, j, i-1, j);
        }
        else if (INU(i,j+1) && INU(i,j-1)) {
          u[i][j] = (u[i][j+1] + u[i][j-1])/2.;
          printf("average u: %d,%d + %d,%d\n", i, j+1, i, j-1);
        }
        else {
        */
#ifdef EXTRA_LINEAR
        extra_velocity((real)i, (real)j + 0.5, u, v, ap, in, 
                       &eu, &ev, nx, ny);
#else
        extra_velocity2((real)i, (real)j + 0.5, u, v, ap, in, 
                        &eu, &ev, nx, ny);
#endif
        u[i][j] = eu;
        /*      } */
      }
      else if (!INU(i,j))
        u[i][j] = UNDEFINED;    
    for (j = 3; j < ny; j++)
      if (!INV(i,j) && BV(i,j)) {
        /*
        if (INV(i+1,j) && INV(i-1,j)) {
          v[i][j] = (v[i+1][j] + v[i-1][j])/2.;
          printf("average v: %d,%d + %d,%d\n", i+1, j, i-1, j);
        }
        else if (INV(i,j+1) && INV(i,j-1)) {
          v[i][j] = (v[i][j+1] + v[i][j-1])/2.;
          printf("average v: %d,%d + %d,%d\n", i, j+1, i, j-1);
        }
        else {
        */
#ifdef EXTRA_LINEAR
          extra_velocity((real)i + 0.5, (real)j, u, v, ap, in, 
                         &eu, &ev, nx, ny);
#else
          extra_velocity2((real)i + 0.5, (real)j, u, v, ap, in, 
                          &eu, &ev, nx, ny);
#endif
          v[i][j] = ev;
          /*    } */
      }
      else if (!INV(i,j))
        v[i][j] = UNDEFINED;
  }
  bc_vector_bound(u, v, nx, ny);
}


void bc_pressure_inside(real2D p, real2D ap, int nx, int ny)
{
  int i, j;
  
  for (i = 2; i < nx; i++)
    for (j = 2; j < ny; j++)
      if (!INP(i,j))
        p[i][j] = 0.0;
}


void bc_pressure(real2D p, real2D ap, real pe, 
                 real R,
                 int nx, int ny)
{
  int i, j;
  
  bc_pressure_inside(p, ap, nx, ny);
  for (i = 2; i < nx; i++) {
#ifdef WALL_TOP
    p[i][ny] = p[i][ny-1];
#else
    p[i][ny] = rp_pressure(pe, R, (real)i + 0.5, (real)ny + 0.5, nx, ny);
#endif
    p[i][1] = p[i][2]; /* axisymmetric */
  }
  for (j = 2; j < ny; j++) {
#ifdef WALL_LEFT
    p[1][j] = p[2][j];
#else
    p[1][j] = rp_pressure(pe, R, 1.5, (real)j + 0.5, nx, ny);
#endif
#ifdef WALL_RIGHT
    p[nx][j] = p[nx-1][j];
#else
    p[nx][j] = rp_pressure(pe, R, (real)nx + 0.5, (real)j + 0.5, nx, ny);
#endif
  }
}


void bc_scalar(real2D scal, int nx, int ny, char sw)
{
  int i, j;
  
  /*  div: gradient equal to zero on y=0,1 planes  */
  switch(sw) {  
    case NULGRAD:
      for (i = 1; i <= nx; i++) {
        scal[i][1] = scal[i][2]; /* axisymmetric */
        scal[i][ny] = scal[i][ny - 1];
      }     
      for (j = 1; j <= ny; j++) {
        scal[1][j] = scal[2][j];
        scal[nx][j] = scal[nx-1][j];
      }
      break;
    case NULGRAD2:
      for (i = 1; i <= nx; i++)
        scal[i][1] = scal[i][2]; /* axisymmetric */
#ifdef WALL_LEFT
      for (j = 1; j <= ny; j++)
        scal[1][j] = scal[2][j];
#endif
#ifdef WALL_RIGHT
      for (j = 1; j <= ny; j++)
        scal[nx][j] = scal[nx-1][j];
#endif
#ifdef WALL_TOP
      for (i = 1; i <= nx; i++)
        scal[i][ny] = scal[i][ny-1];
#endif
      break;
    case NUL:
      /*  residue equal to zero on the walls  */ 
      for (i = 1; i <= nx; i++) {
        scal[i][1] = 0.0; /* axisymmetric */
        scal[i][ny] = 0.0;
      }
      for (j = 1; j <= ny; j++)
        scal[1][j] = scal[nx][j] = 0.0;
      break;
    default:
      printf("illegal switch value for a scalar field\n");
      exit(1);
  }
}


int bc_xcoord(int i, int nx)
{
  /* mirror symmetric boundary conditions */
  if (i <= 1)
    return 3 - i;
  if (i >= nx)
    return 2*nx - i - 1;
  return i;
}


int bc_ycoord(int j, int ny)
{
  /* mirror symmetric boundary conditions */
  if (j <= 1)
    return 3 - j;
  if (j >= ny)
    return 2*ny - j - 1;
  return j;
}

---