00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024 static char const Id[] = "$Id: ipol.c 2962 2010-02-04 19:50:44Z scholz $\n\n";
00025 static char const Td[] = "$Today: " __FILE__ " " __DATE__ " " __TIME__ " $\n\n";
00026
00027 #include "util.h"
00028
00029 #define MAXLINELEN 100000
00030 #define HTCOL 3
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055 int readht(FILE *fd,PetscReal **indata,const int col)
00056 {
00057 int nrows;
00058 PetscReal *xydata;
00059
00060 MagparFunctionLogBegin;
00061
00062 int rank,size;
00063 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
00064 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
00065
00066 if (col<2) {
00067 SETERRQ1(PETSC_ERR_ARG_CORRUPT,"col=%i < 2!\n",col);
00068 }
00069
00070
00071 if (!rank) {
00072
00073 fscanf(fd,"%i\n",&nrows);
00074 PetscPrintf(PETSC_COMM_WORLD,"Reading column %i, %i lines\n",col,nrows);
00075
00076
00077 nrows += 2;
00078
00079 ierr = PetscMalloc(HTCOL*nrows*sizeof(PetscReal),&xydata);CHKERRQ(ierr);
00080
00081
00082 for (int i=1; i<nrows-1; i++) {
00083 int k;
00084
00085 #ifdef PETSC_USE_SINGLE
00086 k=fscanf(fd,"%f",xydata+i*HTCOL+0);
00087 #else
00088 k=fscanf(fd,"%lf",xydata+i*HTCOL+0);
00089 #endif
00090 if (k!=1) {
00091 SETERRQ2(PETSC_ERR_FILE_UNEXPECTED,
00092 "xydata corrupt in row %i col %i\n",
00093 i,1
00094 );
00095 }
00096
00097
00098 for (int j=1; j<col; j++) {
00099 #ifdef PETSC_USE_SINGLE
00100 k=fscanf(fd,"%f",xydata+i*HTCOL+1);
00101 #else
00102 k=fscanf(fd,"%lf",xydata+i*HTCOL+1);
00103 #endif
00104
00105 if (k!=1) {
00106 SETERRQ2(PETSC_ERR_FILE_UNEXPECTED,
00107 "xydata in corrupt in row %i col %i\n",
00108 i,j
00109 );
00110 }
00111 }
00112
00113 char msg[MAXLINELEN];
00114 fgets(msg,MAXLINELEN,fd);
00115 }
00116
00117
00118 xydata[0]=PETSC_MIN;
00119 xydata[1]=0.0;
00120 xydata[2]=0.0;
00121
00122
00123 xydata[HTCOL*(nrows-1)+0]=PETSC_MAX;
00124 xydata[HTCOL*(nrows-1)+1]=0.0;
00125 xydata[HTCOL*(nrows-1)+2]=0.0;
00126
00127
00128 for (int i=0; i<nrows-1; i++) {
00129 xydata[i*HTCOL+2] = (xydata[(i+1)*HTCOL+1]-xydata[i*HTCOL+1])/
00130 (xydata[(i+1)*HTCOL+0]-xydata[i*HTCOL+0]);
00131 }
00132
00133 PetscPrintf(PETSC_COMM_WORLD,
00134 "%14s %14s %14s\n",
00135 "time","ts1","dts1/dt"
00136 );
00137
00138 for (int i=0; i<nrows; i++) {
00139 PetscPrintf(PETSC_COMM_WORLD,
00140 "%14e %14e %14e\n",
00141 xydata[HTCOL*i+0],
00142 xydata[HTCOL*i+1],
00143 xydata[HTCOL*i+2]
00144 );
00145 }
00146 }
00147
00148
00149 ierr = MPI_Bcast(&nrows,1,MPI_INT,0,PETSC_COMM_WORLD);CHKERRQ(ierr);
00150 if (rank) {
00151 ierr = PetscMalloc(HTCOL*nrows*sizeof(PetscReal),&xydata);CHKERRQ(ierr);
00152 }
00153 ierr = MPI_Bcast(xydata,nrows*HTCOL,MPIU_SCALAR,0,PETSC_COMM_WORLD);CHKERRQ(ierr);
00154 *indata=xydata;
00155
00156 MagparFunctionLogReturn(0);
00157 }
00158
00159 PetscReal getscalf(PetscReal xydata[],PetscReal xval)
00160 {
00161 int tint;
00162 PetscReal result;
00163
00164 tint=0;
00165 while (xval > xydata[HTCOL*(tint+1)]) tint++;
00166
00167
00168 result=(xval-xydata[HTCOL*tint+0])*xydata[HTCOL*tint+2]+xydata[HTCOL*tint+1];
00169
00170 return(result);
00171 }
00172
00173 int readxy(FILE *fd,const int colx,const int coly,int *n_rows,PetscReal **indata)
00174 {
00175 int nrows;
00176 PetscReal *xydata;
00177
00178 MagparFunctionLogBegin;
00179
00180 int rank,size;
00181 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
00182 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
00183
00184 if (colx<1) {
00185 SETERRQ1(PETSC_ERR_ARG_CORRUPT,"colx=%i < 1!\n",colx);
00186 }
00187 if (coly<1) {
00188 SETERRQ1(PETSC_ERR_ARG_CORRUPT,"coly=%i < 1!\n",coly);
00189 }
00190
00191
00192 if (!rank) {
00193
00194 fscanf(fd,"%i\n",&nrows);
00195 PetscPrintf(PETSC_COMM_WORLD,"Reading columns %i, %i, and %i lines\n",colx,coly,nrows);
00196
00197 ierr = PetscMalloc(2*nrows*sizeof(PetscReal),&xydata);CHKERRQ(ierr);
00198
00199
00200 for (int i=0;i<nrows;i++) {
00201 PetscReal val,x,y;
00202 x=y=val=PETSC_MIN;
00203
00204
00205 for (int j=1;j<=PetscMax(colx,coly);j++) {
00206
00207
00208 int k;
00209 #ifdef PETSC_USE_SINGLE
00210 k=fscanf(fd,"%f",&val);
00211 #else
00212 k=fscanf(fd,"%lf",&val);
00213 #endif
00214 if (j==colx && k==1) x=val;
00215 if (j==coly && k==1) y=val;
00216
00217 if (x!=PETSC_MIN && y!=PETSC_MIN) {
00218 xydata[2*i+0]=x;
00219 xydata[2*i+1]=y;
00220 break;
00221 }
00222 }
00223
00224
00225 char msg[MAXLINELEN];
00226 fgets(msg,MAXLINELEN,fd);
00227 }
00228
00229 #ifdef DEBUG
00230 PetscPrintf(PETSC_COMM_WORLD,"%14s %14s\n","x","y");
00231
00232 for (int i=0;i<nrows;i++) {
00233 PetscPrintf(PETSC_COMM_WORLD,
00234 "%14e %14e\n",
00235 xydata[2*i+0],
00236 xydata[2*i+1]
00237 );
00238 }
00239 #endif
00240 }
00241
00242
00243 ierr = MPI_Bcast(&nrows,1,MPI_INT,0,PETSC_COMM_WORLD);CHKERRQ(ierr);
00244 if (rank) {
00245 ierr = PetscMalloc(2*nrows*sizeof(PetscReal),&xydata);CHKERRQ(ierr);
00246 }
00247 ierr = MPI_Bcast(xydata,nrows*2,MPIU_SCALAR,0,PETSC_COMM_WORLD);CHKERRQ(ierr);
00248 *indata=xydata;
00249
00250 *n_rows=nrows;
00251
00252 MagparFunctionLogReturn(0);
00253 }
00254