CDIwrite.c 7.25 KB
Newer Older
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
/*
  This file is part of CDO. CDO is a collection of Operators to
  manipulate and analyse Climate model Data.

  Copyright (C) 2003-2012 Uwe Schulzweida, Uwe.Schulzweida@zmaw.de
  See COPYING file for copying and redistribution conditions.

  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; version 2 of the License.

  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.
*/


#include <cdi.h>
#include "cdo.h"
#include "cdo_int.h"
#include "pstream.h"


Uwe Schulzweida's avatar
Uwe Schulzweida committed
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
static
const char *filetypestr(int filetype)
{
  switch ( filetype )
    {
    case FILETYPE_GRB:  return ("GRIB");            break;
    case FILETYPE_GRB2: return ("GRIB2");           break;
    case FILETYPE_NC:   return ("netCDF");          break;
    case FILETYPE_NC2:  return ("netCDF2");         break;
    case FILETYPE_NC4:  return ("netCDF4");         break;
    case FILETYPE_NC4C: return ("netCDF4 classic"); break;
    case FILETYPE_SRV:  return ("SERVICE");         break;
    case FILETYPE_EXT:  return ("EXTRA");           break;
    case FILETYPE_IEG:  return ("IEG");             break;
    default:            return ("");
    }
}

static
const char *datatypestr(int datatype)
{
  static char str[20];

  str[0] = 0;
  sprintf(str, "%d bit packed", datatype);

  if      ( datatype == DATATYPE_PACK   ) return ("P0");
  else if ( datatype > 0 && datatype <= 32 ) return (str);
  else if ( datatype == DATATYPE_CPX32  ) return ("C32");
  else if ( datatype == DATATYPE_CPX64  ) return ("C64");
  else if ( datatype == DATATYPE_FLT32  ) return ("32 bit floats");
  else if ( datatype == DATATYPE_FLT64  ) return ("64 bit floats");
  else if ( datatype == DATATYPE_INT8   ) return ("I8");
  else if ( datatype == DATATYPE_INT16  ) return ("I16");
  else if ( datatype == DATATYPE_INT32  ) return ("I32");
  else if ( datatype == DATATYPE_UINT8  ) return ("U8");
  else if ( datatype == DATATYPE_UINT16 ) return ("U16");
  else if ( datatype == DATATYPE_UINT32 ) return ("U32");
  else                                    return ("");
}

static
off_t filesize(const char *filename)
{
  FILE *fp;
  off_t pos = 0;

  fp = fopen(filename, "r");
  if ( fp == NULL )
    {
      fprintf(stderr, "Open failed on %s\n", filename);
    }
  else
    {
      fseek(fp, 0L, SEEK_END);
      pos = ftello(fp);
    }
  
  return pos;
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
86
87
void *CDIwrite(void *argument)
{
88
  int memtype = MEMTYPE_DOUBLE;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
89
90
  int nvars = 10, nlevs = 0, ntimesteps = 30;
  char *defaultgrid = "global_.2";
Uwe Schulzweida's avatar
Uwe Schulzweida committed
91
92
  int operatorID;
  int streamID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
93
  int tsID, varID, levelID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
94
  int gridsize, i;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
95
  int rval, rstart, rinc;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
96
97
98
  int vlistID;
  int gridID = -1, zaxisID, taxisID;
  int vdate, vtime, julday;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
99
100
  int filetype, datatype;
  off_t fsize;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
101
  unsigned int seed = 1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
102
  const char *gridfile;
103
  char *envstr;
104
  off_t values = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
105
  double file_size, data_size = 0;
106
  double tw, t0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
107
108
  double *levels = NULL;
  double ***vars = NULL;
109
  float *farray = NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
110
111
112
  extern int timer_write;

  srand(seed);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
113
114
115

  cdoInitialize(argument);

116
117
118
119
120
121
122
  envstr = getenv("MEMTYPE");
  if ( envstr )
    {
      if      ( strcmp(envstr, "float")  == 0 ) memtype = MEMTYPE_FLOAT;
      else if ( strcmp(envstr, "double") == 0 ) memtype = MEMTYPE_DOUBLE;
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
123
124
125
  if ( cdoVerbose ) cdoPrint("parameter: <grid, <nlevs, <ntimesteps, <nvars>>>>");
  // operatorInputArg("<grid, <nlevs, <ntimesteps, <nvars>>>>");

Uwe Schulzweida's avatar
Uwe Schulzweida committed
126
  if ( operatorArgc() > 4 ) cdoAbort("Too many arguments!");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
127

Uwe Schulzweida's avatar
Uwe Schulzweida committed
128
129
  gridfile = defaultgrid;
  if ( operatorArgc() >= 1 ) gridfile = operatorArgv()[0];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
130
131
132
  if ( operatorArgc() >= 2 ) nlevs = atol(operatorArgv()[1]);
  if ( operatorArgc() >= 3 ) ntimesteps = atol(operatorArgv()[2]);
  if ( operatorArgc() >= 4 ) nvars = atol(operatorArgv()[3]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
133

Uwe Schulzweida's avatar
Uwe Schulzweida committed
134
  gridID   = cdoDefineGrid(gridfile);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
135
  gridsize = gridInqSize(gridID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
136

Uwe Schulzweida's avatar
Uwe Schulzweida committed
137
  zaxisID  = zaxisCreate(ZAXIS_SURFACE, 1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
138

Uwe Schulzweida's avatar
Uwe Schulzweida committed
139
  if ( cdoVerbose )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
140
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
141
142
143
144
145
146
147
148
149
150
151
152
      cdoPrint("gridsize   : %d", gridInqSize);
      cdoPrint("nlevs      : %d", nlevs);
      cdoPrint("ntimesteps : %d", ntimesteps);
      cdoPrint("nvars      : %d", nvars);
    } 

  if ( nlevs <= 0 ) nlevs = 1;
  if ( ntimesteps <= 0 ) ntimesteps = 1;
  if ( nvars <= 0 ) nvars = 1;

  vars = (double ***) malloc(nvars*sizeof(double **));
  for ( varID = 0; varID < nvars; varID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
153
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
154
155
156
157
158
159
160
      vars[varID] = (double **) malloc(nlevs*sizeof(double *));
      for ( levelID = 0; levelID < nlevs; levelID++ )
	{
	  vars[varID][levelID] = (double *) malloc(gridsize*sizeof(double));
	  for ( i = 0; i < gridsize; ++i )
	    vars[varID][levelID][i] = varID + rand()/(RAND_MAX+1.0);
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
161
162
    }

163
164
  if ( memtype == MEMTYPE_FLOAT ) farray = (float *) malloc(gridsize*sizeof(float));

Uwe Schulzweida's avatar
Uwe Schulzweida committed
165
166
  vlistID = vlistCreate();

Uwe Schulzweida's avatar
Uwe Schulzweida committed
167
  for ( i = 0; i < nvars; ++i )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
168
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
169
170
171
      varID = vlistDefVar(vlistID, gridID, zaxisID, TIME_VARIABLE);
      vlistDefVarCode(vlistID, varID, varID+1);
      //    vlistDefVarName(vlistID, varID, );
Uwe Schulzweida's avatar
Uwe Schulzweida committed
172
173
174
175
176
    }

  taxisID = taxisCreate(TAXIS_RELATIVE);
  vlistDefTaxis(vlistID, taxisID);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
177
  // vlistDefNtsteps(vlistID, 1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
178
179
180
181
182

  streamID = streamOpenWrite(cdoStreamName(0), cdoFiletype());

  streamDefVlist(streamID, vlistID);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
183
184
185
186
  filetype = streamInqFiletype(streamID);
  datatype = vlistInqVarDatatype(vlistID, 0);
	  
  julday = date_to_julday(CALENDAR_PROLEPTIC, 19870101);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
187

188
189
  t0 = timer_val(timer_write);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
190
191
192
193
194
195
196
197
198
199
200
  for ( tsID = 0; tsID < ntimesteps; tsID++ )
    {
      rval  = rstart + rinc*tsID;
      vdate = julday_to_date(CALENDAR_PROLEPTIC, julday + tsID);
      vtime = 0;
      taxisDefVdate(taxisID, vdate);
      taxisDefVtime(taxisID, vtime);
      streamDefTimestep(streamID, tsID);

      for ( varID = 0; varID < nvars; varID++ )
        {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
201
          for ( levelID = 0; levelID < nlevs; levelID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
202
            {
203
	      values += gridsize;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
204
              streamDefRecord(streamID, varID, levelID);
205
206
	      if ( memtype == MEMTYPE_FLOAT )
		{
207
208
		  double *darray = vars[varID][levelID];
		  for ( i = 0; i < gridsize; ++i ) farray[i] = darray[i];
209
210
211
212
213
214
215
216
		  streamWriteRecordFloat(streamID, farray, 0);
		  data_size += gridsize*4;
		}
	      else
		{
		  streamWriteRecord(streamID, vars[varID][levelID], 0);
		  data_size += gridsize*8;
		}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
217
218
            }
        }
219
220
221
222
223
224
225

      if ( cdoVerbose )
	{
	  tw = timer_val(timer_write) - t0;
	  t0 = timer_val(timer_write);
	  cdoPrint("Timestep %d: %.2f seconds", tsID+1, tw);
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
226
227
228
229
    }

  streamClose(streamID);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
230
231
  tw = timer_val(timer_write);

232
  values /= 1000000;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
233
  data_size /= 1024.*1024.*1024.;
234
  if ( memtype == MEMTYPE_FLOAT )
235
    cdoPrint("Wrote %.1f GB of 32 bit floats to %s %s, %.1f milval/s", data_size, datatypestr(datatype), filetypestr(filetype), values/tw);
236
  else
237
    cdoPrint("Wrote %.1f GB of 64 bit floats to %s %s, %.1f milval/s", data_size, datatypestr(datatype), filetypestr(filetype), values/tw);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
238
239
240
241
242
243

  fsize = filesize(cdoStreamName(0));
  file_size = fsize;
  file_size /= 1024.*1024.*1024.;
  cdoPrint("Wrote %.1f GB in %.1f seconds, total %.1f MB/s", file_size, tw, 1024*file_size/tw);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
244
245
  vlistDestroy(vlistID);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
246
247
248
  for ( varID = 0; varID < nvars; varID++ )
    {
      for ( levelID = 0; levelID < nlevs; levelID++ ) free(vars[varID][levelID]);
249
      free(vars[varID]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
250
251
    }
  free(vars);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
252

253
254
  if ( farray ) free(farray);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
255
256
257
258
  cdoFinish();

  return (0);
}