---

fas.c

Go to the documentation of this file.

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

#include "utilf.h"
#include "mymalloc.h"
#include "make_bc.h"
#include "fas.h"

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

#define COARSEST     1
#define COARSERELAX  4

static real2D  *itau, *irhs, *irho, *ip, *ia, *ib, *iei, *itemp, *icc;

static void abfromab(real2D a, real2D b, real2D ei, real2D a2, real2D c2, 
                     int n1, int n2);

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

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

static void slvsml(real2D u, real2D rhs, 
                   real2D a, real2D b, real2D ei, 
                   int nfx, int nfy);

void lop(real2D out, real2D u, real2D a, real2D b, 
         int n1, int n2);

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

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;
  
  return n;
}

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

int initpressure(int n1, int n2, int ng)
{
  int j;
  
  irhs = (real2D *) malloc((ng + 1) * sizeof(real2D));
  irho = (real2D *) malloc((ng + 1) * sizeof(real2D));
  itau = (real2D *) malloc((ng + 1) * sizeof(real2D));
  ip = (real2D *) malloc((ng + 1) * sizeof(real2D));
  ia = (real2D *) malloc((ng + 1) * sizeof(real2D));
  ib = (real2D *) malloc((ng + 1) * sizeof(real2D));
  iei = (real2D *) malloc((ng + 1) * sizeof(real2D));
  itemp = (real2D *) malloc((ng + 1) * sizeof(real2D));
  icc = (real2D *) malloc((ng + 1) * sizeof(real2D));
  
  if (icc EQUALS NULL)
      return MEMERROR;
  
  irhs[ng] = (real2D)mymalloc(sizeof(real), n1, n2);
  itau[ng] = (real2D)mymalloc(sizeof(real), n1, n2); 
  itemp[ng] = (real2D)mymalloc(sizeof(real), n1, n2);

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

  for (j = ng - 1; j >= 1; j--)
    {
      irhs[j] = (real2D)mymalloc(sizeof(real), n1, n2); 
      irho[j] = (real2D)mymalloc(sizeof(real), n1, n2); 
      itau[j] = (real2D)mymalloc(sizeof(real), n1, n2); 
      ip[j] = (real2D)mymalloc(sizeof(real), n1, n2); 
      ia[j] = (real2D)mymalloc(sizeof(real), n1, n2);
      ib[j] = (real2D)mymalloc(sizeof(real), n1, n2); 
      iei[j] = (real2D)mymalloc(sizeof(real), n1, n2); 
      icc[j] = (real2D)mymalloc(sizeof(real), n1, n2);
      itemp[j] = (real2D)mymalloc(sizeof(real), n1, n2);

      if (itemp[j] EQUALS NULL)
        return MEMERROR;

      n1 = n1 / 2 + 1;
      n2 = n2 / 2 + 1;
    }
    
  return 0;
} /* end initpressure() */

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

void mginit(real2D u, real2D p, real2D a, real2D b, real2D ei, real2D cc,
            int n1, int n2, int ng) 
{
  int ngrid, nn1, nn2;
  
  ip[ng] = p;
  irho[ng] = u;
  ia[ng] = a;
  ib[ng] = b;
  iei[ng] = ei;
  icc[ng] = cc;

  nn1 = n1 / 2 + 1;
  nn2 = n2 / 2 + 1;
  ngrid = ng - 1;
  
  bc_scalar(u, n1, n2, NULGRAD);

  rstrct(irho[ngrid], u, nn1, nn2);
  rstrct(icc[ngrid], p, nn1, nn2);  
  abfromab(ia[ngrid], ib[ngrid], iei[ngrid],
           ia[ngrid+1], ib[ngrid+1], nn1, nn2);
  while ((nn1 > 4) && (nn2 > 4))
    {
      nn1 = nn1 / 2 + 1;
      nn2 = nn2 / 2 + 1;
      ngrid = ngrid - 1;
      rstrct(irho[ngrid], irho[ngrid+1], nn1, nn2);
      rstrct(icc[ngrid], icc[ngrid+1], nn1, nn2);
      abfromab(ia[ngrid], ib[ngrid], iei[ngrid], ia[ngrid+1], ib[ngrid+1],
               nn1, nn2);
    }
}

void mgfas(real2D u, real2D p, real2D a, real2D b, real2D ei, 
           int n1, int n2, int ng, int npre, int npost, int ncycle)
{
  int i, j, jj, jcycle, jpre, jpost, ngrid, nfx, nfy, nn1, nn2;

  ip[ng] = p;
  irho[ng] = u;
  ia[ng] = a;
  ib[ng] = b;
  iei[ng] = ei;  
  icc[ng] = cc;

  nn1 = n1 / 2 + 1;
  nn2 = n2 / 2 + 1;
  ngrid = ng - 1;
  
  bc_scalar(u, n1, n2, NULGRAD);

  rstrct(irho[ngrid], u, nn1, nn2);
  rstrct(icc[ngrid], cc, nn1, nn2);
  abfromab(ia[ngrid], ib[ngrid], iei[ngrid],
           ia[ngrid+1], ib[ngrid+1], nn1, nn2);
  while ((nn1 > 4) && (nn2 > 4))
    {
      nn1 = nn1 / 2 + 1;
      nn2 = nn2 / 2 + 1;
      ngrid = ngrid - 1;      
      rstrct(irho[ngrid], irho[ngrid+1], nn1, nn2);
      rstrct(icc[ngrid], icc[ngrid+1], nn1, nn2);
      abfromab(ia[ngrid], ib[ngrid], iei[ngrid], ia[ngrid+1], ib[ngrid+1],
               nn1, nn2);
    }

  ngrid = ng;
  for (j = 2; j <= ngrid; j++)
    {
      nn1 = 2 * nn1 - 2;
      nn2 = 2 * nn2 - 2;
      interp(ip[j], ip[j-1], nn1, nn2);
      copy(irhs[j], irho[j], nn1, nn2);
      for (jcycle = 1; jcycle <= ncycle; jcycle++)
        {
          nfx = nn1;
          nfy = nn2;
          for (jj = j; jj >= 2; jj--)
            {
              for (jpre = 1; jpre <= npre; jpre++)
                {
                  relax(ip[jj], irhs[jj], ia[jj], ib[jj], iei[jj], icc[jj],
                        nfx, nfy);
                  bc_scalar(ip[jj], nfx, nfy, NULGRAD);
                }
              lop(itemp[jj], ip[jj], ia[jj], ib[jj], icc[jj], nfx, nfy);              
              bc_scalar(itemp[jj], nfx, nfy, NULGRAD);
              nfx = nfx / 2 + 1;
              nfy = nfy / 2 + 1;
              rstrct(itemp[jj-1], itemp[jj], nfx, nfy);
              rstrct(ip[jj-1], ip[jj], nfx, nfy);
              lop(itau[jj-1], ip[jj-1], ia[jj-1], ib[jj-1], icc[jj-1], nfx, nfy);
              bc_scalar(itau[jj-1], nfx, nfy, NULGRAD);
              matsub(itau[jj-1], itemp[jj-1], itau[jj-1], nfx, nfy);
              rstrct(irhs[jj-1], irhs[jj], nfx, nfy);
              matadd(irhs[jj-1], itau[jj-1], irhs[jj-1], nfx, nfy);
            }
          
          slvsml(ip[1], irhs[1],
                 ia[1], ib[1], iei[1], icc[1], nfx, nfy);
          
          for (jj = 2; jj <= j; jj++)
            {
              rstrct(itemp[jj-1], ip[jj], nfx, nfy);
              matsub(ip[jj-1], itemp[jj-1], itemp[jj-1], nfx, nfy);
              nfx = 2 * nfx - 2;
              nfy = 2 * nfy - 2;
              interp(itau[jj], itemp[jj-1], nfx, nfy);
              matadd(ip[jj], itau[jj], ip[jj], nfx, nfy);
              for (jpost = 1; jpost <= npost; jpost++)
                {
                  relax(ip[jj], irhs[jj], ia[jj], ib[jj], iei[jj], icc[jj],
                        nfx, nfy);
                  bc_scalar(ip[jj], nfx, nfy, NULGRAD);
                }
            }  
        }
    }
}
#if 0
/* --------------------------------------------------------------------- */

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

  /* solve exactly on the coarsest grid */
  if (n EQUALS COARSEST)
    {
      slvsml(ip[COARSEST], irhs[COARSEST], 
             ia[COARSEST], ib[COARSEST], iei[COARSEST], icc[COARSEST},
             nx, ny);
      return;
    }
  
  /* do prerelaxation */
  for (count = 0; count < npre; count++)
    {
      relax(ip[n], irhs[n], ia[n], ib[n], iei[n], icc[n], nx, ny);
      bc_scalar(ip[n], nx, ny, NULGRAD);
    }

  /* compute tau = LH(Rph) - RLh(ph) */
  lop(itemp[n], ip[n], ia[n], ib[n], icc[n], nx, ny);
  bc_scalar(itemp[n], nx, ny, NULGRAD);
  rstrct(itemp[n-1], itemp[n], ncx, ncy);
  rstrct(ip[n-1], ip[n], ncx, ncy);
  lop(itau[n-1], ip[n-1], ia[n-1], ib[n-1], icc[n-1], ncx, ncy);
  bc_scalar(itau[n-1], ncx, ncy, NULGRAD);
  matsub(itau[n-1], itemp[n-1], itau[n-1], ncx, ncy);
  
  /* compute r.h.s. */
  rstrct(irhs[n-1], irhs[n], ncx, ncy);
  matadd(irhs[n-1], itau[n-1], irhs[n-1], ncx, ncy);

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

  /* compute tau = P(pH-Rph) */
  rstrct(itemp[n-1], ip[n], ncx, ncy);
  matsub(ip[n-1], itemp[n-1], itemp[n-1], ncx, ncy);
  printf("--------------------itemp[%d]--------------------\n", n-1);
  matprint(itemp[n-1], ncx, ncy);
  interp(itau[n], itemp[n-1], nx, ny);
  printf("--------------------itau[%d]------------------\n", n);
  matprint(itau[n], nx, ny);

  /* correct on the fine grid */
  printf("--------------------ip[%d]--------------------\n", n);
  matprint(ip[n], nx, ny);
  matadd(ip[n], itau[n], ip[n], nx, ny);

  /* do post relaxation */
  for (count = 0; count < npost; count++)
    {
      relax(ip[n], irhs[n], ia[n], ib[n], iei[n], icc[n], nx, ny);
      bc_scalar(ip[n], nx, ny, NULGRAD);
    }
}
#endif
/* --------------------------------------------------------------------- */

static void abfromab(real2D a, real2D b, real2D ei, real2D a2, real2D c2, 
                     int n1, int n2)
{
  int i, j;
  
  for (i = 1; i <= n1 - 1; i++)
    for (j = 2; j <= n2 - 1; j++)
      a[i][j] = 0.5 * (a2[2*i-1][2*j-2] + a2[2*i-1][2*j-1]);
  for (i = 2; i <= n1 - 1; i++)
    for (j = 1; j <= n2 - 1; j++)
      b[i][j] = 0.5 * (c2[2*i-2][2*j-1] + c2[2*i-1][2*j-1]);
  for (i = 2; i <= n1 - 1; i++)
    for (j = 2; j <= n2 - 1; j++)
      ei[i][j] = 1.0 / (a[i-1][j] + a[i][j] + b[i][j-1] + b[i][j]);
}

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

static void rstrct(real2D uc, real2D uf, int ncx, int ncy)
{
  int ic, iff, jc, jf;
  
  for (jc = 2; jc <= ncx - 1; jc++)
    {
      jf = 2 * jc - 2;
      for (ic = 2; ic <= ncy - 1; ic++)
        {
          iff = 2 * ic - 2;
          uc[jc][ic] = 0.5 * uf[jf][iff] + 
            0.125 * (uf[jf][iff+1] + uf[jf][iff-1] + 
                     uf[jf+1][iff] + uf[jf-1][iff]);
        }
    }
  bc_scalar(uc, ncx, ncy, NULGRAD);
}

/* --------------------------------------------------------------------- */
#if 0
static void interp(real2D uf, real2D uc, int nfx, int nfy)
{
  int ic, iif, jc, jf, ncx = nfx / 2 + 1, ncy = nfy / 2 + 1;

  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]);

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

static void interp(real2D uf, real2D uc, int nfx, int nfy)
{
  int ic, iff, jc, jf, ncx, ncy;
  
  ncx = nfx / 2 + 1;
  ncy = nfy / 2 + 1;
  for (ic = 2; ic <= ncx - 1; ic++)
    {
      iff = 2 * ic - 2;
      for (jc = 2; jc <= ncy - 1; jc++)
        {
          jf = 2 * jc - 2;
          uf[iff][jf] = 
          uf[iff][jf+1] = 
          uf[iff+1][jf] = 
          uf[iff+1][jf+1] = uc[ic][jc];
        }
    }
  
  bc_scalar(uf, nfx, nfy, NULGRAD);
}

#if 0
static void interp(real2D uf, real2D uc, int nfx, int nfy)
{
  int ic, iff, jc, jf, ncx, ncy;
  
  ncx = nfx / 2 + 1;
  ncy = nfy / 2 + 1;

  for (ic = 2; ic <= ncx - 1; ic++)
    {
      iff = 2 * ic - 2;
      for (jc = 2; jc <= ncy - 1; jc++)
        {
          jf = 2 * jc - 2;
          uf[iff][jf] = uc[ic][jc];
        }
    }
  
  for (jf = 1; jf <= nfy; jf += 2)
    for (iff = 2; iff < nfx; iff += 2)
      uf[iff][jf] = 0.5 * (uf[iff+1][jf] + uf[iff-1][jf]);
  for (jf = 2; jf < nfy; jf += 2)
    for (iff = 1; iff <= nfx; iff++)
      uf[iff][jf] = 0.5 * (uf[iff][jf+1] + uf[iff][jf-1]);

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

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

static void slvsml(real2D u, real2D rhs, real2D a, real2D b, real2D ei, 
                   int nfx, int nfy)
{
  int count;
  
  for (count = 1; count <= COARSERELAX; count++)
    {
      relax(u, rhs, a, b, ei, nfx, nfy);
      bc_scalar(u, nfx, nfy, NULGRAD);
    }
}

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

void relax(real2D u, real2D rhs, 
           real2D a, real2D b, real2D ei, real2D cc,
           int n1, int n2)
{
  real h2, h, res;
  int i, j, isw, jsw;

  h = 1.0 / (n2 - 2);
  h2 = h * h;
  for (jsw = 1; jsw <= 2; jsw++)
    {
      isw = jsw;
      for (i = 2; i <= n1 - 1; i++, isw = 3 - isw)
        for (j = isw + 1; j <= n2 - 1; j += 2)
          if (cc[i][j] == 1.0)
            u[i][j] = (a[i][j] * u[i+1][j] +
                       a[i-1][j] * u[i-1][j] +
                       b[i][j] * u[i][j+1] +
                       b[i][j-1] * u[i][j-1]
                       - h2 * rhs[i][j]) 
              * ei[i][j];
    }
}


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

void lop(real2D out, real2D u, real2D a, real2D b, real2D cc,
         int n1, int n2)
{
  int i, j;
  real hi2;

  hi2 = (real) (n2 - 2) * (n2 - 2);
  for (i = 2; i <= n1 - 1; i++)
    for (j = 2; j <= n2 - 1; j++)
      if (cc[i][j] == 1.0)
        out[i][j] = hi2 * (a[i][j] * (u[i+1][j] - u[i][j]) +
                           a[i-1][j] * (u[i-1][j] - u[i][j]) +
                           b[i][j] * (u[i][j+1] - u[i][j]) +
                           b[i][j-1] * (u[i][j-1] - u[i][j]));
}

---