psn_1d.c

Finite Volume Poisson PDE Solver

/*psn_lib/psn_1d.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    : 1D POISSON SOLVER
    VERSION        : 1.0
************************************************************************/
/*17.01.04*/

#include "psn.h"


int psn_1d_eps_solve_low(
   const int N,
   const int grid,
   const double x0,
   const double x1,
   const double alpha0,
   const double alpha1,
   const double beta0,
   const double beta1,
   const double gamma0,
   const double gamma1,
   const PSN_VECTOR v,
   const PSN_VECTOR e,
   PSN_VECTOR u,
   PSN_VECTOR a,
   PSN_VECTOR b,
   PSN_VECTOR c
)
{
  int errcode;
  double hx, hx2, tm;
  int i;

  if (grid>0)
    hx=(x1-x0)/N;
  else
    hx=(x1-x0)/(N+1);
  hx2=hx*hx;
  // Initialize tridiagonal matrix

  if (e)
    for (i = 1; i <= N; i++) {
      a[i] = (i>1)? 2*e[i-1]/(e[i-1]+e[i]):0;
      c[i] = (i<N)?2*e[i+1]/(e[i]+e[i+1]):0;
      b[i] = -a[i]-c[i];
      u[i]=v[i]*hx2/e[i];
    }
  else
    for (i = 1; i <= N; i++) {
      a[i] = (i>1)? 1:0;
      c[i] = (i<N)?1:0;
      b[i] = -a[i]-c[i];
      u[i]=v[i]*hx2;
    }

  //Boundary conditions
  if (grid>0) { //staggered grid
    b[1]-=(2.0+4.*beta0 / (alpha0 * hx - 2*beta0));
    b[N]-=(2.0-4.*beta1 / (alpha1 * hx + 2*beta1));
    u[1] -= 2.0*gamma0 * hx / (alpha0 * hx - 2*beta0);
    u[N] -= 2.0*gamma1 * hx / (alpha1 * hx + 2*beta1);
  } else {  //standard grid
    b[1] -=(1+ beta0 / (alpha0 * hx - beta0));
    b[N] -=(1- beta1 / (alpha1 * hx + beta1));
    u[1] -= gamma0 * hx / (alpha0 * hx - beta0);
    u[N] -= gamma1 * hx / (alpha1 * hx + beta1);
  }
  // Invert tridiagonal matrix equation

  errcode=psn_tridiag_solve_low(N,a,b,c,u);

  // Calculate i=0 and i=N+1 values
  tm= (grid>0) ? 2.:1.;
  u[0] = (gamma0 * hx - tm*beta0 * u[1]) /(alpha0 * hx - tm*beta0);
  u[N+1] = (gamma1 * hx + tm*beta1 * u[N]) /(alpha1 * hx + tm*beta1);
  return errcode;
}

int psn_1d_struct_init(psn_1d_struct *pst, int N)
{
   pst->N=N;
   pst->v=psn_vector_create(N+2);
   pst->e=psn_vector_create(N+2);
   pst->u=psn_vector_create(N+2);
   pst->a=psn_vector_create(N+2);
   pst->b=psn_vector_create(N+2);
   pst->c=psn_vector_create(N+2);
   return 0;
}

int psn_1d_struct_free(psn_1d_struct *pst)
{
   pst->N=0;
   psn_vector_free(pst->v);
   psn_vector_free(pst->e);
   psn_vector_free(pst->u);
   psn_vector_free(pst->a);
   psn_vector_free(pst->b);
   psn_vector_free(pst->c);
   return 0;
}

int psn_1d_solver(psn_1d_struct *pst)
{
   if (!pst->N) return 1;
   return psn_1d_eps_solve_low(
      pst->N,
      pst->grid,
      pst->x[0],
      pst->x[1],
      pst->alpha[0],
      pst->alpha[1],
      pst->beta[0],
      pst->beta[1],
      pst->gamma[0],
      pst->gamma[1],
      pst->v,
      pst->e,
      pst->u,
      pst->a,
      pst->b,
      pst->c
   );
}