psn_2d.c

Finite Volume Poisson PDE Solver

/*psn_lib/psn_2d.c*/
/***********************************************************************
    Finite Volume Poisson PDE Solver: C-Library & Matlab Toolbox
    Implements numerical solution of Poisson PDE
    in 2D  Cartesian and Cylindrical coordinates

    Copyright (C) 2004 Igor Kaufman
    Copyright (C) 2004 Lancaster University, UK

    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA


    Author's email: [e-mail]

    APPLICATION    : 2D POISSON SOLVER
    VERSION        : 1.0
************************************************************************/
/*17.01.04*/

#include "psn.h"

#define set_boundval(alpha,beta,gamma,h,u) (gamma*h-beta*(u))/(alpha*h-beta)
#define b_bound(alpha,beta,gamma,h)        (alpha*h)/(alpha*h-beta)
#define w_bound(alpha,beta,gamma,h)        (gamma*h)/(alpha*h-beta);
#define over_relax(u_old, u_new, w)        (1-w)*(u_old)+w*(u_new)
#define W_OPT(Nx, Ny)                       2-6/(max(Nx,Ny))

void set_bounds(psn_2d_struct *psd,  double hx,  double hy )
{

    double alpha, beta, gamma, h, tmp;
    size_t col, row, i, row_b, col_b;
    size_t Nx=psd->Nx;
    size_t Ny=psd->Ny;
    PSN_MATRIX u=psd->u;
    PSN_MATRIX v=psd->v;

    h=hy/(psd->grid+1);
    for (col=1;col<=Nx;col++) {
       for (i=LOW;i<=UP;i++){
          if (i==UP || psd->crd==0) {
            alpha=psd->alphaY[i];
            beta=psd->betaY[i];
            switch (i) {
              case LOW:
                row=1;
                row_b=0;
                break;
              case UP:
                row=Ny;
                row_b=Ny+1;
                beta=-beta;
                break;
            }
            gamma=psn_matrix_getval(v, row_b,col);
            tmp=set_boundval(alpha,beta,gamma,h,psn_matrix_getval(u,row,col));
          } else {
            row_b=0;
            tmp=psn_matrix_getval(u,1,col);
          }
          psn_matrix_setval(u,row_b,col,tmp);
       }
    }

    h=hx/(psd->grid+1);
    for (row=1;row<=Ny;row++) {
       for (i=LOW;i<=UP;i++){
          alpha=psd->alphaX[i];
          beta=psd->betaX[i];
          switch (i) {
            case LOW:
              col=1;
              col_b=0;
              break;
            case UP:
              col=Nx;
              col_b=Nx+1;
              beta=-beta;
              break;
          }
          gamma=psn_matrix_getval(v,row,col_b);
          tmp=set_boundval(alpha,beta,gamma,h,psn_matrix_getval(u, row,col));
          psn_matrix_setval(u,row,col_b,tmp);
       }
    }
    tmp=0.5*(psn_matrix_getval(u,0,1)+psn_matrix_getval(u,1,0));
    psn_matrix_setval(u,0,0,tmp);
    tmp=0.5*(psn_matrix_getval(u,0,Nx)+psn_matrix_getval(u,1,Nx+1));
    psn_matrix_setval(u,0,Nx+1,tmp);
    tmp=0.5*(psn_matrix_getval(u,Ny+1,1)+psn_matrix_getval(u,Ny,0));
    psn_matrix_setval(u,Ny+1,0,tmp);
    tmp=0.5*(psn_matrix_getval(u,Ny,Nx+1)+psn_matrix_getval(u,Ny+1,Nx));
    psn_matrix_setval(u,Ny+1,Nx+1,tmp);

}

void set_diff2_coeff(psn_2d_struct *psd,
            double hx,
            double hy,
            size_t row,
            size_t col,

            double a[],
            double b[],
            double c[],
            double w[]
    )
{
   double alpha, beta, gamma, ee;
   PSN_MATRIX e=psd->e;
   PSN_MATRIX v=psd->v;
   size_t dim;


   //Set interior coefficients
   if  (e==NULL)  { //Vacuum
     a[0]= (col>1)? 1.:0;
     c[0]= (col<psd->Nx)? 1.:0;
     a[1]=(row>1)? 1.:0;
     c[1]=(row<psd->Ny)? 1.:0;
   } else {  //Matter
     ee=psn_matrix_getval(e, row,col);
     a[0]=(col>1)? 2.*psn_matrix_getval(e,row,col-1)/(psn_matrix_getval(e,row,col-1)+ee):0;
     c[0]=(col<psd->Nx)? 2.*psn_matrix_getval(e,row,col+1)/(psn_matrix_getval(e,row,col+1)+ee):0;
     a[1]=(row>1) ? 2.*psn_matrix_getval(e,row-1,col)/(psn_matrix_getval(e,row-1,col)+ee):0;
     c[1]=(row<psd->Ny)? 2.*psn_matrix_getval(e,row+1,col)/(psn_matrix_getval(e,row+1,col)+ee):0;
   }
   if (psd->crd==1)  {  //cylindrical coordinates
     double t=1.0/(2*row-1);
     a[1]*=(2.*(row-1))*t;
     c[1]*=(2.*row)*t;
   }
   for (dim=0;dim<2;dim++) {
     b[dim]=-a[dim]-c[dim];
     w[dim]=0;
   }


   if (col==1 || col==psd->Nx) { //Add boundary conditions for X=0 and X=Lx
     double h=hx/(psd->grid+1);
     if (col==1) {
       alpha=psd->alphaX[0];
       beta=psd->betaX[0];
       gamma=psn_matrix_getval(v,row,0);
     } else {
       alpha=psd->alphaX[1];
       beta=-psd->betaX[1];
       gamma=psn_matrix_getval(v,row,psd->Nx+1);
     }
     b[0]-=b_bound(alpha,beta,gamma,h);
     w[0]=w_bound(alpha,beta,gamma,h);
   }

   if (row==1 || row==psd->Ny) {  //Add boundary conditions for Y=0 and Y=Ly
     double h=hy/(psd->grid+1);
     if (row==1) {
       if (psd->crd==0) {   //Don't do for cylindrical where R==0
         alpha=psd->alphaY[0];
         beta=psd->betaY[0];
         gamma=psn_matrix_getval(v,0,col);
       }
     } else {
       alpha=psd->alphaY[1];
       beta=-psd->betaY[1];
       gamma=psn_matrix_getval(v,psd->Ny+1, col);
     }
     if (psd->crd==0 || row==psd->Ny) {
       b[1]-=b_bound(alpha,beta,gamma,h);
       w[1]=w_bound(alpha,beta,gamma,h);
     }
   }

}

int psn_2d_struct_init(psn_2d_struct *psd, const size_t Nx, const size_t Ny)
{
   psd->Nx=Nx;
   psd->Ny=Ny;
   psd->v=psn_matrix_create(Ny,Nx);
   psd->e=psn_matrix_create(Ny,Nx);
   psd->u=psn_matrix_create(Ny,Nx);
   return 0;
}

int psn_2d_struct_free(psn_2d_struct *psd)
{
   psn_matrix_free(psd->v,psd->Ny);
   psn_matrix_free(psd->e,psd->Ny);
   psn_matrix_free(psd->u,psd->Ny);
   psd->Nx=0;
   psd->Ny=0;
   return 0;
}

int psn_2d_solver_LSOR(psn_2d_struct *psd)
{

   int res;
   int stop;
   size_t step, row, col;
   size_t Nx=psd->Nx;
   size_t Ny=psd->Ny;
   PSN_MATRIX u=psd->u, v=psd->v, e=psd->e;

   size_t Nm=(Nx>Ny)?Nx:Ny;
   double hx, hy, W, h2[2];
   double err;

   PSN_VECTOR um=psn_vector_create(Nm+2);
   PSN_VECTOR wm=psn_vector_create(Nm+2);
   PSN_VECTOR a=psn_vector_create(Nm+2);
   PSN_VECTOR b=psn_vector_create(Nm+2);
   PSN_VECTOR c=psn_vector_create(Nm+2);

   //Set grid
   if (psd->grid) {
     hx=psd->Lx/(Nx);
     hy=psd->Ly/(Ny);
   } else {
     hx=psd->Lx/(Nx+1);
     hy=psd->Ly/(Ny+1);
   }
   h2[0]=1/(hx*hx);
   h2[1]=1/(hy*hy);


   W=psd->W;
   if (W==0) W=2.-6./max(Nx,Ny);

   //Main iteration cycle
   for (step=0, stop=0;step<psd->max_step && !stop;step++) {
     err=0;
     for (row=1;row<=Ny;row++) {
        psn_matrix_get_rowcol(Nx,row,smRow,u,um);  //save old row for U
        //Define coeff
        for (col=1;col<=Nx;col++) {
          double aa[2],bb[2], cc[2], ww[2];
          set_diff2_coeff(psd,hx,hy,row,col, aa,bb,cc,ww);
          a[col]=h2[0]*aa[0];
          b[col]=h2[0]*bb[0]+h2[1]*bb[1];
          c[col]=h2[0]*cc[0];
          wm[col]=psn_matrix_getval(v, row,col);
          if (e) wm[col]/=psn_matrix_getval(e, row,col);
          wm[col]-=h2[0]*ww[0];
          wm[col]-=h2[1]*(ww[1]+aa[1]*psn_matrix_getval(u, row-1,col)+cc[1]*psn_matrix_getval(u,row+1,col));
        }
        //and solve for current row
        psn_tridiag_solve_low(Nx,a,b,c,wm);
        //OverRelaxation
        if (W!=1)
          for (col=1;col<=Nx;col++)
             wm[col]=over_relax(um[col],wm[col],W);

        err=max(err,(norm_vector(Nx,wm,um)));
        psn_matrix_set_rowcol(Nx,row,smRow,u,wm);
     }

     if (psd->method==metADI)
     for (col=1;col<=Nx;col++) {
        psn_matrix_get_rowcol(Nx,col,smCol,u,um);  //save old row for U
        //Define coeff
        for (row=1;row<=Ny;row++) {
          double aa[2],bb[2], cc[2], ww[2];
          set_diff2_coeff(psd,hx,hy,row,col,aa,bb,cc,ww);
          a[row]=h2[1]*aa[1];
          b[row]=h2[1]*bb[1]+h2[0]*bb[0];
          c[row]=h2[1]*cc[1];
          wm[row]=psn_matrix_getval(v, row,col);
          if (e) wm[row]/=psn_matrix_getval(e, row,col);
          wm[row]-=h2[1]*ww[1];
          wm[row]-=h2[0]*(ww[0]+aa[0]*psn_matrix_getval(u, row,col-1)+cc[0]*psn_matrix_getval(u,row,col+1));
        }
        //and solve for current row
        psn_tridiag_solve_low(Ny,a,b,c,wm);
        //OverRelaxation
        if (W!=1)
          for (row=1;row<=Ny;row++)
             wm[row]=over_relax(um[row],wm[row],W);

        //err=max(err,(norm_vector(Ny,wm,um)));
        psn_matrix_set_rowcol(Ny,col,smCol,u,wm);
     }

     if (psd->step_err) psd->step_err[step]=err;
     stop=(err<psd->tol);
   }

   set_bounds(psd,hx,hy);
   psd->error=err;
   psd->num_step=step;
   res=(psd->error<psd->tol)? 0:3;
   psd->result=res;

   psn_vector_free(um);
   psn_vector_free(wm);
   psn_vector_free(a);
   psn_vector_free(b);
   psn_vector_free(c);

   return res;
}

int psn_2d_solver_PSOR(psn_2d_struct *psd)
{

   int res;
   int stop;
   size_t step, row, col;
   size_t Nx=psd->Nx;
   size_t Ny=psd->Ny;
   PSN_MATRIX u=psd->u, v=psd->v, e=psd->e;

   double hx, hy, W;
   double kappa, kappa2, hv;
   double err;

   //Set grid
   if (psd->grid) {
     hx=psd->Lx/(Nx);
     hy=psd->Ly/(Ny);
   } else {
     hx=psd->Lx/(Nx+1);
     hy=psd->Ly/(Ny+1);
   }
   kappa=hx/hy;
   kappa2=kappa*kappa;
   hv=hx*hx;


   W=psd->W;
   if (W==0) W=2.-6./max(Nx,Ny);

   //Main iteration cycle
   for (step=0, stop=0;step<psd->max_step && !stop;step++) {
      double aa[2], bb[2], cc[2], ww[2], rhs, b01;
      double r, uval, delta;

      err=0;
      for (row=1;row<=Ny;row++) {
        for (col=1;col<=Nx;col++) {
          set_diff2_coeff(psd,hx,hy,row,col,aa,bb,cc,ww);
          b01=bb[0]+bb[1]*kappa2;
          rhs=psn_matrix_getval(v,row,col)*hv;
          if (e) rhs/=psn_matrix_getval(e, row,col);
          uval=psn_matrix_getval(u,row,col);
          ww[0]+=(aa[0]*psn_matrix_getval(u,row,col-1)+cc[0]*psn_matrix_getval(u,row,col+1));
          ww[1]+=(aa[1]*psn_matrix_getval(u,row-1,col)+cc[1]*psn_matrix_getval(u,row+1,col));
          r=ww[0]+ww[1]*kappa2+b01*uval-rhs;
          delta=-W*r/b01;
          err=max(err, fabs(delta));
          psn_matrix_setval(u,row,col,(uval+delta));
        }
      }

      if (psd->step_err) psd->step_err[step]=err;
      stop=(err<psd->tol);
   }
   set_bounds(psd,hx,hy);
   psd->error=err;
   psd->num_step=step;
   res=(psd->error<psd->tol)? 0:3;
   psd->result=res;

   return res;
}

int psn_2d_solver(psn_2d_struct *psd)
{
   switch (psd->method) {
      case metLSOR:
      case metADI:
        return psn_2d_solver_LSOR(psd);
      case metPSOR:
        return psn_2d_solver_PSOR(psd);
      default:
        return psn_2d_solver_PSOR(psd);
   }
}