ipol.c

Go to the documentation of this file.

/*
    This file is part of magpar.

    Copyright (C) 2002-2010 Werner Scholz

    www:   http://www.magpar.net/
    email: magpar(at)magpar.net

    magpar 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.

    magpar 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 magpar; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

static char const Id[] = "$Id: ipol.c 2962 2010-02-04 19:50:44Z scholz $\n\n";
static char const Td[] = "$Today: " __FILE__ "  " __DATE__ "  "  __TIME__  " $\n\n";

#include "util.h"

#define MAXLINELEN 100000
#define HTCOL 3

/*
input file structure:
-----------
n
x_0  y_0
x_1  y_1
x_2  y_2
x_3  y_3
...
x_n  y_n
-----------

xydata structure:
-----------------------------
x_0  y_0  (y_1-y_0)/(x_1-x_0)
x_1  y_1  (y_2-y_1)/(x_2-x_1)
x_2  y_2  (y_3-y_2)/(x_3-x_2)
x_3  y_3  (y_4-y_3)/(x_4-x_3)
...
x_n  y_n  0
-----------------------------
*/

int readht(FILE *fd,PetscReal **indata,const int col)
{
  int            nrows;
  PetscReal      *xydata;

  MagparFunctionLogBegin;

  int rank,size;
  ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);

  if (col<2) {
    SETERRQ1(PETSC_ERR_ARG_CORRUPT,"col=%i < 2!\n",col);
  }

  /* only the first processor reads the file */
  if (!rank) {
    /* get number of lines */
    fscanf(fd,"%i\n",&nrows);
    PetscPrintf(PETSC_COMM_WORLD,"Reading column %i, %i lines\n",col,nrows);

    /* account for additional first and last entry */
    nrows += 2;

    ierr = PetscMalloc(HTCOL*nrows*sizeof(PetscReal),&xydata);CHKERRQ(ierr);

    /* get xydata from file (leave room for first ) */
    for (int i=1; i<nrows-1; i++) {
      int k;
      /* read x from first column */
#ifdef PETSC_USE_SINGLE
      k=fscanf(fd,"%f",xydata+i*HTCOL+0);
#else
      k=fscanf(fd,"%lf",xydata+i*HTCOL+0);
#endif
      if (k!=1) {
        SETERRQ2(PETSC_ERR_FILE_UNEXPECTED,
          "xydata corrupt in row %i col %i\n",
          i,1
        );
      }

      /* read y from column col */
      for (int j=1; j<col; j++) {
#ifdef PETSC_USE_SINGLE
        k=fscanf(fd,"%f",xydata+i*HTCOL+1);
#else
        k=fscanf(fd,"%lf",xydata+i*HTCOL+1);
#endif

        if (k!=1) {
          SETERRQ2(PETSC_ERR_FILE_UNEXPECTED,
            "xydata in corrupt in row %i col %i\n",
            i,j
          );
        }
      }
      /* read remainder of this line */
      char msg[MAXLINELEN];
      fgets(msg,MAXLINELEN,fd);
    }

    /* set first xydata set */
    xydata[0]=PETSC_MIN;
    xydata[1]=0.0;
    xydata[2]=0.0;

    /* set last xydata set */
    xydata[HTCOL*(nrows-1)+0]=PETSC_MAX;
    xydata[HTCOL*(nrows-1)+1]=0.0;
    xydata[HTCOL*(nrows-1)+2]=0.0;

    /* calculate dH1/dt and dH2/dt */
    for (int i=0; i<nrows-1; i++) {
      xydata[i*HTCOL+2] = (xydata[(i+1)*HTCOL+1]-xydata[i*HTCOL+1])/
                          (xydata[(i+1)*HTCOL+0]-xydata[i*HTCOL+0]);
    }

    PetscPrintf(PETSC_COMM_WORLD,
      "%14s %14s %14s\n",
      "time","ts1","dts1/dt"
    );

    for (int i=0; i<nrows; i++) {
      PetscPrintf(PETSC_COMM_WORLD,
        "%14e %14e %14e\n",
        xydata[HTCOL*i+0],
        xydata[HTCOL*i+1],
        xydata[HTCOL*i+2]
      );
    }
  }

  /* broadcast xydata to all processors */
  ierr = MPI_Bcast(&nrows,1,MPI_INT,0,PETSC_COMM_WORLD);CHKERRQ(ierr);
  if (rank) {
    ierr = PetscMalloc(HTCOL*nrows*sizeof(PetscReal),&xydata);CHKERRQ(ierr);
  }
  ierr = MPI_Bcast(xydata,nrows*HTCOL,MPIU_SCALAR,0,PETSC_COMM_WORLD);CHKERRQ(ierr);
  *indata=xydata;

  MagparFunctionLogReturn(0);
}

PetscReal getscalf(PetscReal xydata[],PetscReal xval)
{
  int       tint;
  PetscReal result;

  tint=0;
  while (xval > xydata[HTCOL*(tint+1)]) tint++;

  /* calculate scaling factor */
  result=(xval-xydata[HTCOL*tint+0])*xydata[HTCOL*tint+2]+xydata[HTCOL*tint+1];

  return(result);
}

int readxy(FILE *fd,const int colx,const int coly,int *n_rows,PetscReal **indata)
{
  int            nrows;
  PetscReal      *xydata;

  MagparFunctionLogBegin;

  int rank,size;
  ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);

  if (colx<1) {
    SETERRQ1(PETSC_ERR_ARG_CORRUPT,"colx=%i < 1!\n",colx);
  }
  if (coly<1) {
    SETERRQ1(PETSC_ERR_ARG_CORRUPT,"coly=%i < 1!\n",coly);
  }

  /* only the first processor reads the file */
  if (!rank) {
    /* get number of lines */
    fscanf(fd,"%i\n",&nrows);
    PetscPrintf(PETSC_COMM_WORLD,"Reading columns %i, %i, and %i lines\n",colx,coly,nrows);

    ierr = PetscMalloc(2*nrows*sizeof(PetscReal),&xydata);CHKERRQ(ierr);

    /* loop over all rows */
    for (int i=0;i<nrows;i++) {
      PetscReal val,x,y;
      x=y=val=PETSC_MIN;

      /* loop over columns */
      for (int j=1;j<=PetscMax(colx,coly);j++) {

        /* read value */
        int k;
#ifdef PETSC_USE_SINGLE
        k=fscanf(fd,"%f",&val);
#else
        k=fscanf(fd,"%lf",&val);
#endif
        if (j==colx && k==1) x=val;
        if (j==coly && k==1) y=val;

        if (x!=PETSC_MIN && y!=PETSC_MIN) {
          xydata[2*i+0]=x;
          xydata[2*i+1]=y;
          break;
        }
      }

      /* read remainder of this line */
      char msg[MAXLINELEN];
      fgets(msg,MAXLINELEN,fd);
    }

#ifdef DEBUG
    PetscPrintf(PETSC_COMM_WORLD,"%14s %14s\n","x","y");

    for (int i=0;i<nrows;i++) {
      PetscPrintf(PETSC_COMM_WORLD,
        "%14e %14e\n",
        xydata[2*i+0],
        xydata[2*i+1]
      );
    }
#endif
  }

  /* broadcast xydata to all processors */
  ierr = MPI_Bcast(&nrows,1,MPI_INT,0,PETSC_COMM_WORLD);CHKERRQ(ierr);
  if (rank) {
    ierr = PetscMalloc(2*nrows*sizeof(PetscReal),&xydata);CHKERRQ(ierr);
  }
  ierr = MPI_Bcast(xydata,nrows*2,MPIU_SCALAR,0,PETSC_COMM_WORLD);CHKERRQ(ierr);
  *indata=xydata;

  *n_rows=nrows;

  MagparFunctionLogReturn(0);
}