---

vmg.c

Go to the documentation of this file.

/* --------------------------------------------------------------------- */
/* 
   SCCS information: %W% %G%
*/
/* --------------------------------------------------------------------- */

#include "utilf.h"
#include "mymalloc.h"
#include "markers.h"
#include "surfaces.h"
#include "extra.h"
#include "make_bc.h"
#include "momentum.h"
#include "mg.h"
#include "inout.h"

/* --------------------------------------------------------------------- */

#define COARSEST     3
#define COARSERELAX  20

static real2D *irhs, *ip, *iap, *ierr;
static real2D *ic0, *ic1, *ic2, *ic3, *ic4;
static real2D sxp, syp, syu, au, cu, sxv, av, cv;
static real2D rz1dr, rz2dr, rr1dz, rr2dz;
static real2D w1, w2, dummy;
static int Ngrid;

static void rstrct(real2D uc, real2D uf, int ncx, int ncy, real scaling);

static void interp(real2D uf, real2D uc, int nfx, int nfy);

static void addint(real2D uf, real2D uc, real2D res, real2D ap,
                   int nfx, int nfy);

static void slvsml(real2D u, real2D rhs, real2D cc, 
                   real2D c0, real2D c1, real2D c2, real2D c3, real2D c4,
                   int nfx, int nfy);

/* --------------------------------------------------------------------- */

int computeng(int n1, int n2)
{
  int nmin, n, nn, j;

  nmin = MIN(n1, n2);
  
  j = nmin - 2;  n = 1;
  while (j > 2) 
  {  j /= 2;   n++;  }
  
  nn = 1;  
  for (j = 1; j <= n; j++) 
    nn *= 2;       
  if((nn + 2) != nmin)
     return NGERROR;

  if ((n1 - 2) % (nmin - 2) != 0)
     return NXERROR;
  if ((n2 - 2) % (nmin - 2) != 0)
     return NYERROR;

  Ngrid = n;
  return n;
}

/* --------------------------------------------------------------------- */

int initpressure(int n1, int n2, int ng)
{
  int j;
  int err = 0;

  err |= ((irhs = (real2D *) malloc((ng + 1) * sizeof(real2D))) == NULL);
  err |= ((ip = (real2D *) malloc((ng + 1) * sizeof(real2D))) == NULL);
  err |= ((iap = (real2D *) malloc((ng + 1) * sizeof(real2D))) == NULL);
  err |= ((ic0 = (real2D *) malloc((ng + 1) * sizeof(real2D))) == NULL);
  err |= ((ic1 = (real2D *) malloc((ng + 1) * sizeof(real2D))) == NULL);
  err |= ((ic2 = (real2D *) malloc((ng + 1) * sizeof(real2D))) == NULL);
  err |= ((ic3 = (real2D *) malloc((ng + 1) * sizeof(real2D))) == NULL);
  err |= ((ic4 = (real2D *) malloc((ng + 1) * sizeof(real2D))) == NULL);
  err |= ((ierr = (real2D *) malloc((ng + 1) * sizeof(real2D))) == NULL);

  if ( err ) return MEMERROR;

  err |= ((sxp = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL); 
  err |= ((syp = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL); 

  err |= ((syu = (real2D)mymalloc(sizeof(real), n1 + 1, n2)) == NULL); 
  err |= ((au = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL); 
  err |= ((cu = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL); 

  err |= ((sxv = (real2D)mymalloc(sizeof(real), n1, n2 + 1)) == NULL);
  err |= ((av = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL); 
  err |= ((cv = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL); 

  err |= ((rz1dr = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL); 
  err |= ((rz2dr = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL); 
  err |= ((rr1dz = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL); 
  err |= ((rr2dz = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL); 

  err |= ((w1 = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL); 
  err |= ((w2 = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL);
  err |= ((dummy = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL); 

  err |= ((ierr[ng] = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL); 
  err |= ((irhs[ng] = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL);

  n1 = n1 / 2 + 1;
  n2 = n2 / 2 + 1;

  for (j = ng - 1; j >= 1; j--) {
    err |= ((irhs[j] = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL);
    err |= ((ip[j] = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL);
    err |= ((iap[j] = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL);
    err |= ((ic0[j] = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL);
    err |= ((ic1[j] = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL);
    err |= ((ic2[j] = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL);
    err |= ((ic3[j] = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL);
    err |= ((ic4[j] = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL);
    err |= ((ierr[j] = (real2D)mymalloc(sizeof(real), n1, n2)) == NULL);
    
    if ( err ) return MEMERROR;
    
    n1 = n1 / 2 + 1;
    n2 = n2 / 2 + 1;
  }

  return 0;
} /* end initpressure() */

/* --------------------------------------------------------------------- */

int mgsolve(real2D u, real2D p, real2D ap,
            real2D c0, real2D c1, real2D c2, real2D c3, real2D c4,
            interface in, real lredis, real maxdiv,
            int n1, int n2, int ng) 
{
  int ngrid, nn1, nn2, i, j;
  interface intemp;
  interface_cut cutp, cutu, cutv;
  real scale = 1.0;
  int ncycle = 0;

  resid(irhs[ng], p, u, ap, c0, c1, c2, c3, c4, n1, n2);
  bc_scalar(irhs[ng], n1, n2, NUL);
  if (rdivmax(irhs[ng], n1, n2, &i, &j) <= maxdiv)
    return 0;

  ip[ng] = p;
  iap[ng] = ap;
  ic0[ng] = c0; ic1[ng] = c1; ic2[ng] = c2; ic3[ng] = c3; ic4[ng] = c4;
  interface_init(&intemp);
  intemp.bc = in.bc;
  cutp.list = cutu.list = cutv.list = NULL;

  nn1 = n1;
  nn2 = n2;
  ngrid = ng;
  while (ngrid > COARSEST) {
    nn1 = nn1 / 2 + 1;
    nn2 = nn2 / 2 + 1;
    ngrid--;
    scale *= 2.0;
    lredis *= 2.0;
    /* rescale the interface */
    interface_redis_scale(&intemp, in, lredis, 1./scale);
    interface_splines(intemp, nn1, nn2);
    /* compute volume fraction */
    rstrct(iap[ngrid], iap[ngrid+1], nn1, nn2, 1.0);

    cut_interface(intemp, &cutp, 0.0, 0.0);
    cut_interface(intemp, &cutu, -0.5, 0.0);
    cut_interface(intemp, &cutv, 0.0, -0.5);

    interface_surfaces(intemp, cutp, sxp, syp, iap[ngrid], dummy,
                       w1, w2, nn1, nn2);
    for (i = 2; i <= nn1; i++)
      for (j = 2; j <= nn2; j++) {
        au[i][j] = 0.5*(iap[ngrid][i][j] + iap[ngrid][i-1][j]);
        av[i][j] = 0.5*(iap[ngrid][i][j] + iap[ngrid][i][j-1]);
      } 
    interface_surfaces(intemp, cutu, dummy, syu, au, cu,
                       w1, w2, nn1, nn2);
    interface_surfaces(intemp, cutv, sxv, dummy, av, cv,
                       w1, w2, nn1, nn2);
    sx_fill(sxp, iap[ngrid], 0.0, 0.0, intemp, nn1, nn2);
    sy_fill(syp, iap[ngrid], 0.0, 0.0, intemp, nn1, nn2);
    sy_fill(syu, au, -0.5, 0.0, intemp, nn1, nn2);
    for (j = 2; j <= nn2; j++)
      syu[nn1+1][j] = (real)j - 1.5;
    sx_fill(sxv, av, 0.0, -0.5, intemp, nn1, nn2);
    for (i = 2; i <= nn1; i++)
      sxv[i][nn2+1] = (real)nn2 - 1.5;

    interface_rzdr(intemp, cutp, -0.5, rz1dr, nn1, nn2);
    interface_rzdr(intemp, cutp, 0.5, rz2dr, nn1, nn2);
    interface_rrdz(intemp, cutp, -0.5, rr1dz, nn1, nn2);
    interface_rrdz(intemp, cutp, 0.5, rr2dz, nn1, nn2);
    
    coli(sxp, syp, syu, sxv, cu, cv, iap[ngrid], av, 
         rz1dr, rz2dr, rr1dz, rr2dz,
         ic0[ngrid], ic1[ngrid], ic2[ngrid], ic3[ngrid], ic4[ngrid], 
         nn1, nn2, scale);
  }
  interface_free(intemp);
  if (cutp.list) free(cutp.list);
  if (cutu.list) free(cutu.list);
  if (cutv.list) free(cutv.list);

  copy(w1, p, n1, n2);
  while (rdivmax(irhs[ng], n1, n2, &i, &j) > maxdiv && ncycle < NCYCLEMAX) {
    fill0(ip[ng], n1, n2);
    mgvcycle(ng, n1, n2, 1, 1);
    printf("ncycle: %d rdivmax: %g\n", 
           ncycle + 1, rdivmax(irhs[ng], n1, n2, &i, &j));
    add(w1, p, n1, n2);
    ncycle++;
  }
  copy(p, w1, n1, n2);  

  return ncycle;
} /* end mginit() */

/* --------------------------------------------------------------------- */

void mgvcycle(int n, int nx, int ny, int npre, int npost)
{
  int count;
  int ncx = nx/2 + 1, ncy = ny/2 + 1;
char s[256];

  /* solve exactly on the coarsest grid */
  if (n == COARSEST) {
    slvsml(ip[COARSEST], irhs[COARSEST], iap[COARSEST], 
           ic0[COARSEST], ic1[COARSEST], ic2[COARSEST], ic3[COARSEST], 
           ic4[COARSEST],
           nx, ny);
    return;
  }
  
  /* do prerelaxation */
  for (count = 0; count < npre; count++) {
    relax(ip[n], irhs[n], iap[n], 
          ic0[n], ic1[n], ic2[n], ic3[n], ic4[n], nx, ny);
    bc_scalar(ip[n], nx, ny, NULGRAD2);
  }

  /* compute defect */
  resid(irhs[n], ip[n], irhs[n], iap[n],
        ic0[n], ic1[n], ic2[n], ic3[n], ic4[n], 
        nx, ny);
  bc_scalar(irhs[n], nx, ny, NUL);

  /* compute r.h.s. on the coarser grid using defect */
  rstrct(irhs[n-1], irhs[n], ncx, ncy, - 4.0);

  /* use zero as initial guess on the coarser grid */
  fill0(ip[n-1], ncx, ncy);

  /* solve on the coarser grid */
  mgvcycle(n - 1, ncx, ncy, npre + 1, npost + 1);

  /* interpolate the error from coarse to fine grid */
  interp(ierr[n], ip[n-1], nx, ny);

  /* new residual res(new) = res(old)-A*error */
  resid(irhs[n], ierr[n], irhs[n], iap[n], 
        ic0[n], ic1[n], ic2[n], ic3[n], ic4[n],
        nx, ny);
  bc_scalar(irhs[n], nx, ny, NUL);

  /* correct e(ngrid) with error = interp(e(ngrid-1)) */
  add(ip[n], ierr[n], nx, ny);

  /* do post relaxation with initial guess = 0 */
  fill0(ierr[n], nx, ny);
  for (count = 0; count < npost; count++) {
    relax(ierr[n], irhs[n], iap[n], 
          ic0[n], ic1[n], ic2[n], ic3[n], ic4[n], nx, ny);
    bc_scalar(ierr[n], nx, ny, NULGRAD);
  }

  /* add correction to correction to solution */
  add(ip[n], ierr[n], nx, ny);

  /* compute residual */
  resid(irhs[n], ierr[n], irhs[n], iap[n],
        ic0[n], ic1[n], ic2[n], ic3[n], ic4[n],
        nx, ny);
  bc_scalar(irhs[n], nx, ny, NUL);
/*
  sprintf(s, "ip%d", n);
  datbarray(ip[n], nx, ny, s);
  sprintf(s, "ipc%d", n);
  sprintf(s, "iap%d", n);
  datbarray(iap[n], nx, ny, s);
*/
}

/* --------------------------------------------------------------------- */

static void rstrct(real2D uc, real2D uf, int ncx, int ncy, real scaling)
{
  int ic, iff, jc, jf;
  
  scaling *= 0.5;
  for (jc = 2; jc <= ncx - 1; jc++)
    {
      jf = 2 * jc - 2;
      for (ic = 2; ic <= ncy - 1; ic++)
        {
          iff = 2 * ic - 2;
#if 0
          uc[jc][ic] = (uf[jf][iff] + 
                        0.25 * (uf[jf][iff+1] + uf[jf][iff-1] + 
                                uf[jf+1][iff] + uf[jf-1][iff])) * scaling;
#else
          uc[jc][ic] = .5 * scaling * (uf[jf][iff] + uf[jf + 1][iff] 
                                       + uf[jf][iff+1] + uf[jf + 1][iff+1]);
#endif
        }
    }
  bc_scalar(uc, ncx, ncy, NULGRAD);
}

/* --------------------------------------------------------------------- */

static void interp(real2D uf, real2D uc, int nfx, int nfy)
{
  int ic, iif, jc, jf, ncx, ncy;

  ncx = nfx / 2 + 1;
  ncy = nfy / 2 + 1;
#if 0
  for (jc = 1, jf = 1; jc < ncy; jc++, jf += 2)
    for (ic = 1; ic < ncx; ic++)
      uf[2*ic-1][jf] = uc[ic][jc];
  for (jf = 1; jf < nfy; jf += 2)
    for (iif = 2; iif < nfx; iif += 2)
      uf[iif][jf] = 0.5 * (uf[iif+1][jf] + uf[iif-1][jf]);
  for (jf = 2; jf < nfy; jf += 2)
    for (iif = 1; iif < nfx; iif++)
      uf[iif][jf] = 0.5 * (uf[iif][jf+1] + uf[iif][jf-1]);
#else
  for (ic = 2; ic <= ncx - 1; ic++) {
    iif = 2 * ic - 2;
    for (jc = 2; jc <= ncy - 1; jc++) {
      jf = 2 * jc - 2;
      uf[iif][jf] =
      uf[iif][jf+1] =
      uf[iif + 1][jf] =
      uf[iif + 1][jf+1] = uc[ic][jc];
    }
  }
#endif

  bc_scalar(uf, nfx, nfy, NULGRAD);
}

/* --------------------------------------------------------------------- */

static void addint(real2D uf, real2D uc, real2D res, real2D ap, 
                   int nfx, int nfy)
{
  int i, j;
  
  interp(res, uc, nfx, nfy); 
  for (i = 2; i < nfx; i++)
    for (j = 2; j < nfy; j++)
      if (INP(i,j))
        uf[i][j] += res[i][j];
}

/* --------------------------------------------------------------------- */

static void slvsml(real2D u, real2D rhs, real2D ap, 
                   real2D c0, real2D c1, real2D c2, real2D c3, real2D c4,
                   int nfx, int nfy)
{
  int count;
  
  for (count = 1; count <= COARSERELAX; count++) {
    relax(u, rhs, ap, c0, c1, c2, c3, c4, nfx, nfy);
    bc_scalar(u, nfx, nfy, NULGRAD2);
  }
}

/* --------------------------------------------------------------------- */

void resid(real2D res, real2D p, real2D rhs, real2D ap, 
           real2D c0, real2D c1, real2D c2, real2D c3, real2D c4,
           int nx, int ny)
{
  int i, j;
  
  for (i = 2; i < nx; i++)
    for (j = 2; j < ny; j++)
      if (INP(i,j))
        res[i][j] = c1[i][j]*p[i+1][j] + c2[i][j]*p[i][j+1] 
          + c3[i][j]*p[i-1][j] + c4[i][j]*p[i][j-1] - c0[i][j]*p[i][j]
          - rhs[i][j];
      else
        res[i][j] = 0.0;
}

/* --------------------------------------------------------------------- */

void relax(real2D p, real2D rhs, real2D ap,
           real2D c0, real2D c1, real2D c2, real2D c3, real2D c4,
           int nx, int ny)
{
  int i, j, isw, jsw;

  for (jsw = 1; jsw <= 2; jsw++) {
    isw = jsw;
    for (i = 2; i < nx; i++, isw = 3 - isw)
      for (j = isw + 1; j < ny; j += 2)
        if (INP(i,j))
          p[i][j] = (c1[i][j]*p[i+1][j] + c2[i][j]*p[i][j+1] +
                     c3[i][j]*p[i-1][j] + c4[i][j]*p[i][j-1] 
                     - rhs[i][j])/c0[i][j];
  }
}

---