/[escript]/trunk/escript/src/DataFactory.cpp
ViewVC logotype

Contents of /trunk/escript/src/DataFactory.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1748 - (show annotations)
Wed Sep 3 06:10:39 2008 UTC (10 years, 11 months ago) by ksteube
File size: 11811 byte(s)
MPI parallelism for Data().dump and load.  Use multiple NetCDF
files, one file per MPI process

1
2 /* $Id$ */
3
4 /*******************************************************
5 *
6 * Copyright 2003-2007 by ACceSS MNRF
7 * Copyright 2007 by University of Queensland
8 *
9 * http://esscc.uq.edu.au
10 * Primary Business: Queensland, Australia
11 * Licensed under the Open Software License version 3.0
12 * http://www.opensource.org/licenses/osl-3.0.php
13 *
14 *******************************************************/
15
16 #include "DataFactory.h"
17 #include "esysUtils/esys_malloc.h"
18
19 #include <boost/python/extract.hpp>
20 #include <iostream>
21 #include <exception>
22 #ifdef USE_NETCDF
23 #include <netcdfcpp.h>
24 #endif
25 #ifdef PASO_MPI
26 #include <mpi.h>
27 #endif
28
29 using namespace boost::python;
30
31 namespace escript {
32
33 Data
34 Scalar(double value,
35 const FunctionSpace& what,
36 bool expanded)
37 {
38 //
39 // an empty shape is a scalar
40 DataArrayView::ShapeType shape;
41 return Data(value,shape,what,expanded);
42 }
43
44 Data
45 Vector(double value,
46 const FunctionSpace& what,
47 bool expanded)
48 {
49 DataArrayView::ShapeType shape(1,what.getDomain().getDim());
50 return Data(value,shape,what,expanded);
51 }
52
53 Data
54 Tensor(double value,
55 const FunctionSpace& what,
56 bool expanded)
57 {
58 DataArrayView::ShapeType shape(2,what.getDomain().getDim());
59 return Data(value,shape,what,expanded);
60 }
61
62 Data
63 Tensor3(double value,
64 const FunctionSpace& what,
65 bool expanded)
66 {
67 DataArrayView::ShapeType shape(3,what.getDomain().getDim());
68 return Data(value,shape,what,expanded);
69 }
70
71 Data
72 Tensor4(double value,
73 const FunctionSpace& what,
74 bool expanded)
75 {
76 DataArrayView::ShapeType shape(4,what.getDomain().getDim());
77 return Data(value,shape,what,expanded);
78 }
79
80 Data
81 load(const std::string fileName,
82 const AbstractDomain& domain)
83 {
84 #ifdef USE_NETCDF
85 NcAtt *type_att, *rank_att, *function_space_type_att;
86 // netCDF error handler
87 NcError err(NcError::silent_nonfatal);
88 int mpi_iam=0, mpi_num=1;
89 // Create the file.
90 #ifdef PASO_MPI
91 MPI_Comm_rank(MPI_COMM_WORLD, &mpi_iam);
92 MPI_Comm_size(MPI_COMM_WORLD, &mpi_num);
93 #endif
94 char *newFileName = Escript_MPI_appendRankToFileName(fileName.c_str(), mpi_num, mpi_iam);
95 NcFile dataFile(newFileName, NcFile::ReadOnly);
96 if (!dataFile.is_valid())
97 throw DataException("Error - load:: opening of netCDF file for input failed.");
98 /* recover function space */
99 if (! (function_space_type_att=dataFile.get_att("function_space_type")) )
100 throw DataException("Error - load:: cannot recover function_space_type attribute from escript netCDF file.");
101 int function_space_type = function_space_type_att->as_int(0);
102 delete function_space_type_att;
103 /* test if function space id is valid and create function space instance */
104 if (! domain.isValidFunctionSpaceType(function_space_type) )
105 throw DataException("Error - load:: function space type code in netCDF file is invalid for given domain.");
106 FunctionSpace function_space=FunctionSpace(domain, function_space_type);
107 /* recover rank */
108 if (! (rank_att=dataFile.get_att("rank")) )
109 throw DataException("Error - load:: cannot recover rank attribute from escript netCDF file.");
110 int rank = rank_att->as_int(0);
111 delete rank_att;
112 if (rank<0 || rank>DataArrayView::maxRank)
113 throw DataException("Error - load:: rank in escript netCDF file is greater than maximum rank.");
114 /* recover type attribute */
115 int type=-1;
116 if ((type_att=dataFile.get_att("type")) ) {
117 char* type_str = type_att->as_string(0);
118 if (strncmp(type_str, "constant", strlen("constant")) == 0 ) {
119 type =0;
120 } else if (strncmp(type_str, "tagged", strlen("tagged")) == 0 ) {
121 type =1;
122 } else if (strncmp(type_str, "expanded", strlen("expanded")) == 0 ) {
123 type =2;
124 }
125 esysUtils::free(type_str);
126 } else {
127 if (! (type_att=dataFile.get_att("type_id")) )
128 throw DataException("Error - load:: cannot recover type attribute from escript netCDF file.");
129 type=type_att->as_int(0);
130 }
131 delete type_att;
132
133 /* recover dimension */
134 int ndims=dataFile.num_dims();
135 int ntags =0 , num_samples =0 , num_data_points_per_sample =0, d=0, len_data_point=1;
136 NcDim *d_dim, *tags_dim, *num_samples_dim, *num_data_points_per_sample_dim;
137 /* recover shape */
138 DataArrayView::ShapeType shape;
139 long dims[DataArrayView::maxRank+2];
140 if (rank>0) {
141 if (! (d_dim=dataFile.get_dim("d0")) )
142 throw DataException("Error - load:: unable to recover d0 from netCDF file.");
143 d=d_dim->size();
144 shape.push_back(d);
145 dims[0]=d;
146 len_data_point*=d;
147 }
148 if (rank>1) {
149 if (! (d_dim=dataFile.get_dim("d1")) )
150 throw DataException("Error - load:: unable to recover d1 from netCDF file.");
151 d=d_dim->size();
152 shape.push_back(d);
153 dims[1]=d;
154 len_data_point*=d;
155 }
156 if (rank>2) {
157 if (! (d_dim=dataFile.get_dim("d2")) )
158 throw DataException("Error - load:: unable to recover d2 from netCDF file.");
159 d=d_dim->size();
160 shape.push_back(d);
161 dims[2]=d;
162 len_data_point*=d;
163 }
164 if (rank>3) {
165 if (! (d_dim=dataFile.get_dim("d3")) )
166 throw DataException("Error - load:: unable to recover d3 from netCDF file.");
167 d=d_dim->size();
168 shape.push_back(d);
169 dims[3]=d;
170 len_data_point*=d;
171 }
172 /* recover stuff */
173 Data out;
174 NcVar *var, *ids_var, *tags_var;
175 if (type == 0) {
176 /* constant data */
177 if ( ! ( (ndims == rank && rank >0) || ( ndims ==1 && rank == 0 ) ) )
178 throw DataException("Error - load:: illegal number of dimensions for constant data in netCDF file.");
179 if (rank == 0) {
180 if (! (d_dim=dataFile.get_dim("l")) )
181 throw DataException("Error - load:: unable to recover d0 for scalar constant data in netCDF file.");
182 int d0 = d_dim->size();
183 if (! d0 == 1)
184 throw DataException("Error - load:: d0 is expected to be one for scalar constant data in netCDF file.");
185 dims[0]=1;
186 }
187 out=Data(0,shape,function_space);
188 if (!(var = dataFile.get_var("data")))
189 throw DataException("Error - load:: unable to find data in netCDF file.");
190 if (! var->get(&(out.getDataPoint(0,0).getData()[0]), dims) )
191 throw DataException("Error - load:: unable to recover data from netCDF file.");
192 } else if (type == 1) {
193 /* tagged data */
194 if ( ! (ndims == rank + 1) )
195 throw DataException("Error - load:: illegal number of dimensions for tagged data in netCDF file.");
196 if (! (tags_dim=dataFile.get_dim("num_tags")) )
197 throw DataException("Error - load:: unable to recover number of tags from netCDF file.");
198 ntags=tags_dim->size();
199 dims[rank]=ntags;
200 int *tags = (int *) esysUtils::malloc(ntags*sizeof(int));
201 if (! ( tags_var = dataFile.get_var("tags")) )
202 {
203 esysUtils::free(tags);
204 throw DataException("Error - load:: unable to find tags in netCDF file.");
205 }
206 if (! tags_var->get(tags, ntags) )
207 {
208 esysUtils::free(tags);
209 throw DataException("Error - load:: unable to recover tags from netCDF file.");
210 }
211
212 DataVector data(len_data_point * ntags, 0., len_data_point * ntags);
213 if (!(var = dataFile.get_var("data")))
214 {
215 esysUtils::free(tags);
216 throw DataException("Error - load:: unable to find data in netCDF file.");
217 }
218 if (! var->get(&(data[0]), dims) )
219 {
220 esysUtils::free(tags);
221 throw DataException("Error - load:: unable to recover data from netCDF file.");
222 }
223 out=Data(DataArrayView(data,shape,0),function_space);
224 for (int t=1; t<ntags; ++t) {
225 out.setTaggedValueFromCPP(tags[t],DataArrayView(data,shape,t*len_data_point));
226 }
227 esysUtils::free(tags);
228 } else if (type == 2) {
229 /* expanded data */
230 if ( ! (ndims == rank + 2) )
231 throw DataException("Error - load:: illegal number of dimensions for exanded data in netCDF file.");
232 if ( ! (num_samples_dim = dataFile.get_dim("num_samples") ) )
233 throw DataException("Error - load:: unable to recover number of samples from netCDF file.");
234 num_samples = num_samples_dim->size();
235 if ( ! (num_data_points_per_sample_dim = dataFile.get_dim("num_data_points_per_sample") ) )
236 throw DataException("Error - load:: unable to recover number of data points per sample from netCDF file.");
237 num_data_points_per_sample=num_data_points_per_sample_dim->size();
238 // check shape:
239 if ( ! (num_samples == function_space.getNumSamples() && num_data_points_per_sample == function_space.getNumDataPointsPerSample()) )
240 throw DataException("Error - load:: data sample layout of file does not match data layout of function space.");
241 // get ids
242 if (! ( ids_var = dataFile.get_var("id")) )
243 throw DataException("Error - load:: unable to find reference ids in netCDF file.");
244 const int* ids_p=function_space.borrowSampleReferenceIDs();
245 int *ids_of_nc = (int *)esysUtils::malloc(num_samples*sizeof(int));
246 if (! ids_var->get(ids_of_nc, (long) num_samples) )
247 {
248 esysUtils::free(ids_of_nc);
249 throw DataException("Error - load:: unable to recover ids from netCDF file.");
250 }
251 // check order:
252 int failed=-1, local_failed=-1, i;
253 #pragma omp parallel private(local_failed)
254 {
255 local_failed=-1;
256 #pragma omp for private(i) schedule(static)
257 for (i=0;i < num_samples; ++i) {
258 if (ids_of_nc[i]!=ids_p[i]) local_failed=i;
259 }
260 #pragma omp critical
261 if (local_failed>=0) failed = local_failed;
262 }
263 /* if (failed>=0)
264 {
265 esysUtils::free(ids_of_nc);
266 throw DataException("Error - load:: data ordering in netCDF file does not match ordering of FunctionSpace.");
267 } */
268 // get the data:
269 dims[rank]=num_data_points_per_sample;
270 dims[rank+1]=num_samples;
271 out=Data(0,shape,function_space,true);
272 if (!(var = dataFile.get_var("data")))
273 {
274 esysUtils::free(ids_of_nc);
275 throw DataException("Error - load:: unable to find data in netCDF file.");
276 }
277 if (! var->get(&(out.getDataPoint(0,0).getData()[0]), dims) )
278 {
279 esysUtils::free(ids_of_nc);
280 throw DataException("Error - load:: unable to recover data from netCDF file.");
281 }
282 if (failed>=0) {
283 try {
284 std::cout << "Information - load: start reordering data from netCDF file " << fileName << std::endl;
285 out.borrowData()->reorderByReferenceIDs(ids_of_nc);
286 }
287 catch (std::exception&) {
288 esysUtils::free(ids_of_nc);
289 throw DataException("Error - load:: unable to reorder data in netCDF file.");
290 }
291 }
292 } else {
293 throw DataException("Error - load:: unknown escript data type in netCDF file.");
294 }
295 return out;
296 #else
297 throw DataException("Error - load:: is not compiled with netCFD. Please contact your insatllation manager.");
298 #endif
299 }
300
301 bool
302 loadConfigured()
303 {
304 #ifdef USE_NETCDF
305 return true;
306 #else
307 return false;
308 #endif
309 }
310
311 Data
312 convertToData(const boost::python::object& value,
313 const FunctionSpace& what)
314 {
315 // first we try to extract a Data object from value
316 extract<Data> value_data(value);
317 if (value_data.check()) {
318 Data extracted_data=value_data();
319 if (extracted_data.isEmpty()) {
320 return extracted_data;
321 } else {
322 return Data(extracted_data,what);
323 }
324 } else {
325 return Data(value,what);
326 }
327 }
328
329 } // end of namespace

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26