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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1748 - (hide annotations)
Wed Sep 3 06:10:39 2008 UTC (12 years, 10 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 jgs 82
2 ksteube 1312 /* $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 jgs 474 #include "DataFactory.h"
17 phornby 1628 #include "esysUtils/esys_malloc.h"
18 jgs 82
19     #include <boost/python/extract.hpp>
20 gross 950 #include <iostream>
21 gross 1487 #include <exception>
22 gross 1023 #ifdef USE_NETCDF
23 ksteube 1312 #include <netcdfcpp.h>
24 gross 1023 #endif
25 ksteube 1748 #ifdef PASO_MPI
26     #include <mpi.h>
27     #endif
28 jgs 82
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 gross 950 Data
81     load(const std::string fileName,
82     const AbstractDomain& domain)
83     {
84 gross 1023 #ifdef USE_NETCDF
85 gross 950 NcAtt *type_att, *rank_att, *function_space_type_att;
86     // netCDF error handler
87 gross 1141 NcError err(NcError::silent_nonfatal);
88 ksteube 1748 int mpi_iam=0, mpi_num=1;
89 gross 950 // Create the file.
90 ksteube 1748 #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 gross 950 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 gross 1141 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 phornby 1628 esysUtils::free(type_str);
126 gross 1141 } 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 gross 950 }
131     delete type_att;
132 gross 1141
133 gross 950 /* recover dimension */
134     int ndims=dataFile.num_dims();
135 gross 983 int ntags =0 , num_samples =0 , num_data_points_per_sample =0, d=0, len_data_point=1;
136 gross 967 NcDim *d_dim, *tags_dim, *num_samples_dim, *num_data_points_per_sample_dim;
137 gross 950 /* 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 gross 983 len_data_point*=d;
147 gross 950 }
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 gross 983 len_data_point*=d;
155 gross 950 }
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 gross 983 len_data_point*=d;
163 gross 950 }
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 gross 983 len_data_point*=d;
171 gross 950 }
172     /* recover stuff */
173     Data out;
174 gross 983 NcVar *var, *ids_var, *tags_var;
175 gross 950 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 phornby 1020 throw DataException("Error - load:: illegal number of dimensions for tagged data in netCDF file.");
196 gross 983 if (! (tags_dim=dataFile.get_dim("num_tags")) )
197 phornby 1020 throw DataException("Error - load:: unable to recover number of tags from netCDF file.");
198 gross 950 ntags=tags_dim->size();
199 gross 983 dims[rank]=ntags;
200 phornby 1628 int *tags = (int *) esysUtils::malloc(ntags*sizeof(int));
201 gross 983 if (! ( tags_var = dataFile.get_var("tags")) )
202 phornby 1020 {
203 phornby 1628 esysUtils::free(tags);
204 gross 983 throw DataException("Error - load:: unable to find tags in netCDF file.");
205 phornby 1020 }
206 gross 983 if (! tags_var->get(tags, ntags) )
207 phornby 1020 {
208 phornby 1628 esysUtils::free(tags);
209 phornby 1020 throw DataException("Error - load:: unable to recover tags from netCDF file.");
210     }
211 gross 983
212 gross 1141 DataVector data(len_data_point * ntags, 0., len_data_point * ntags);
213 gross 983 if (!(var = dataFile.get_var("data")))
214 phornby 1020 {
215 phornby 1628 esysUtils::free(tags);
216 phornby 1020 throw DataException("Error - load:: unable to find data in netCDF file.");
217     }
218 gross 983 if (! var->get(&(data[0]), dims) )
219 phornby 1020 {
220 phornby 1628 esysUtils::free(tags);
221 phornby 1020 throw DataException("Error - load:: unable to recover data from netCDF file.");
222     }
223 gross 983 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 phornby 1628 esysUtils::free(tags);
228 gross 950 } 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 gross 967 if ( ! (num_samples_dim = dataFile.get_dim("num_samples") ) )
233 gross 950 throw DataException("Error - load:: unable to recover number of samples from netCDF file.");
234 gross 967 num_samples = num_samples_dim->size();
235     if ( ! (num_data_points_per_sample_dim = dataFile.get_dim("num_data_points_per_sample") ) )
236 gross 950 throw DataException("Error - load:: unable to recover number of data points per sample from netCDF file.");
237 gross 967 num_data_points_per_sample=num_data_points_per_sample_dim->size();
238 gross 982 // check shape:
239 gross 967 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 gross 982 // 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 phornby 1628 int *ids_of_nc = (int *)esysUtils::malloc(num_samples*sizeof(int));
246 gross 982 if (! ids_var->get(ids_of_nc, (long) num_samples) )
247 phornby 1020 {
248 phornby 1628 esysUtils::free(ids_of_nc);
249 phornby 1020 throw DataException("Error - load:: unable to recover ids from netCDF file.");
250     }
251 gross 982 // check order:
252     int failed=-1, local_failed=-1, i;
253     #pragma omp parallel private(local_failed)
254     {
255 gross 985 local_failed=-1;
256 gross 982 #pragma omp for private(i) schedule(static)
257 gross 985 for (i=0;i < num_samples; ++i) {
258 gross 982 if (ids_of_nc[i]!=ids_p[i]) local_failed=i;
259 gross 985 }
260 gross 982 #pragma omp critical
261     if (local_failed>=0) failed = local_failed;
262     }
263 gross 1487 /* if (failed>=0)
264 phornby 1020 {
265 phornby 1628 esysUtils::free(ids_of_nc);
266 phornby 1020 throw DataException("Error - load:: data ordering in netCDF file does not match ordering of FunctionSpace.");
267 gross 1487 } */
268 gross 982 // get the data:
269 gross 967 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 phornby 1020 {
274 phornby 1628 esysUtils::free(ids_of_nc);
275 phornby 1020 throw DataException("Error - load:: unable to find data in netCDF file.");
276     }
277 gross 967 if (! var->get(&(out.getDataPoint(0,0).getData()[0]), dims) )
278 phornby 1020 {
279 phornby 1628 esysUtils::free(ids_of_nc);
280 phornby 1020 throw DataException("Error - load:: unable to recover data from netCDF file.");
281     }
282 gross 1487 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 trankine 1631 catch (std::exception&) {
288 phornby 1628 esysUtils::free(ids_of_nc);
289 gross 1487 throw DataException("Error - load:: unable to reorder data in netCDF file.");
290     }
291     }
292 gross 950 } else {
293     throw DataException("Error - load:: unknown escript data type in netCDF file.");
294     }
295     return out;
296 gross 1023 #else
297     throw DataException("Error - load:: is not compiled with netCFD. Please contact your insatllation manager.");
298     #endif
299     }
300 gross 950
301 gross 1023 bool
302     loadConfigured()
303     {
304     #ifdef USE_NETCDF
305     return true;
306     #else
307     return false;
308     #endif
309 gross 950 }
310    
311 jgs 82 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