Chaste  Release::3.4
CompareHdf5ResultsFiles.cpp
1 /*
2 
3 Copyright (c) 2005-2016, University of Oxford.
4 All rights reserved.
5 
6 University of Oxford means the Chancellor, Masters and Scholars of the
7 University of Oxford, having an administrative office at Wellington
8 Square, Oxford OX1 2JD, UK.
9 
10 This file is part of Chaste.
11 
12 Redistribution and use in source and binary forms, with or without
13 modification, are permitted provided that the following conditions are met:
14  * Redistributions of source code must retain the above copyright notice,
15  this list of conditions and the following disclaimer.
16  * Redistributions in binary form must reproduce the above copyright notice,
17  this list of conditions and the following disclaimer in the documentation
18  and/or other materials provided with the distribution.
19  * Neither the name of the University of Oxford nor the names of its
20  contributors may be used to endorse or promote products derived from this
21  software without specific prior written permission.
22 
23 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29 GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30 HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32 OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33 
34 */
35 
36 #include "CompareHdf5ResultsFiles.hpp"
37 
38 #include <iostream>
39 #include <vector>
40 #include "petscvec.h"
41 #include "Hdf5DataReader.hpp"
42 #include "DistributedVectorFactory.hpp"
43 
44 bool CompareFilesViaHdf5DataReader(std::string pathname1, std::string filename1, bool makeAbsolute1,
45  std::string pathname2, std::string filename2, bool makeAbsolute2,
46  double tol, std::string datasetName)
47 {
48  Hdf5DataReader reader1(pathname1, filename1, makeAbsolute1, datasetName);
49  Hdf5DataReader reader2(pathname2, filename2, makeAbsolute2, datasetName);
50 
51  unsigned number_nodes1 = reader1.GetNumberOfRows();
52  unsigned number_nodes2 = reader2.GetNumberOfRows();
53  if (number_nodes1 != number_nodes2)
54  {
55  std::cout << "Number of nodes " << number_nodes1 << " and " << number_nodes2 << " don't match\n";
56  return false;
57  }
58 
59  // Check the variable names and units
60  std::vector<std::string> variable_names1 = reader1.GetVariableNames();
61  std::vector<std::string> variable_names2 = reader2.GetVariableNames();
62  std::string unlimited_variable1 = reader1.GetUnlimitedDimensionName();
63  std::string unlimited_variable2 = reader2.GetUnlimitedDimensionName();
64  std::string unlimited_variable_unit1 = reader1.GetUnlimitedDimensionUnit();
65  std::string unlimited_variable_unit2 = reader2.GetUnlimitedDimensionUnit();
66 
67  unsigned num_vars = variable_names1.size();
68  if (num_vars != variable_names2.size())
69  {
70  std::cout << "Number of variables " << variable_names1.size()
71  << " and " << variable_names2.size() << " don't match\n";
72  return false;
73  }
74  if (unlimited_variable1 != unlimited_variable2)
75  {
76  std::cout << "Unlimited variable names " << unlimited_variable1 << " and "
77  << unlimited_variable2 << " don't match\n";
78  return false;
79  }
80  if (unlimited_variable_unit1 != unlimited_variable_unit2)
81  {
82  std::cout << "Unlimited variable units " << unlimited_variable_unit1 << " and "
83  << unlimited_variable_unit2 << " don't match\n";
84  return false;
85  }
86  for (unsigned var=0; var<num_vars; var++)
87  {
88  std::string var_name = variable_names1[var];
89  if (var_name != variable_names2[var])
90  {
91  std::cout << "Variable names " << var_name << " and "
92  << variable_names2[var] << " don't match\n";
93  return false;
94  }
95  if (reader1.GetUnit(var_name) != reader2.GetUnit(var_name))
96  {
97  std::cout << "Units names " << reader1.GetUnit(var_name)
98  << " and " << reader2.GetUnit(var_name) << " don't match\n";
99  return false;
100  }
101  }
102 
103  // Check the timestep vectors
104  std::vector<double> times1 = reader1.GetUnlimitedDimensionValues();
105  std::vector<double> times2 = reader2.GetUnlimitedDimensionValues();
106 
107  if (times1.size() != times2.size())
108  {
109  std::cout << "Number of time steps " << times1.size()
110  << " and " << times2.size() << " don't match\n";
111  return false;
112  }
113 
114  for (unsigned timestep=0; timestep<times1.size(); timestep++)
115  {
117  if (fabs(times1[timestep]-times2[timestep]) > 1e-8)
118  {
119  std::cout << "Time steps " << times1[timestep]
120  << " and " << times2[timestep] << " don't match\n";
121  return false;
122  }
123  }
124 
125  bool is_complete1 = reader1.IsDataComplete();
126  bool is_complete2 = reader2.IsDataComplete();
127 
128  if (is_complete1 != is_complete2)
129  {
130  std::cout<<"One of the readers has incomplete data and the other doesn't\n";
131  return false;
132  }
133 
134  if (is_complete1)
135  {
136  DistributedVectorFactory factory(number_nodes1);
137 
138  Vec data1 = factory.CreateVec();
139  Vec data2 = factory.CreateVec();
140 
141  for (unsigned timestep=0; timestep<times1.size(); timestep++)
142  {
143  for (unsigned var=0; var<num_vars; var++)
144  {
145  reader1.GetVariableOverNodes(data1, variable_names1[var], timestep);
146  reader2.GetVariableOverNodes(data2, variable_names2[var], timestep);
147 
148 #if (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 2) //PETSc 2.2
149  double minus_one = -1.0;
150  VecAXPY(&minus_one, data2, data1);
151 #else
152  //[note: VecAXPY(y,a,x) computes y = ax+y]
153  VecAXPY(data1, -1.0, data2);
154 #endif
155 
156  PetscReal difference_norm;
157  VecNorm(data1, NORM_2, &difference_norm);
158 
159  if (difference_norm > tol)
160  {
161  std::cout << "Time " << times1[timestep] << ": vectors differ in NORM_2 by " << difference_norm << std::endl;
162  return false;
163  }
164  }
165  }
166  PetscTools::Destroy(data1);
167  PetscTools::Destroy(data2);
168  }
169  else
170  {
171  // Incomplete data
172 
173  // Check the index vectors
174  std::vector<unsigned> indices1 = reader1.GetIncompleteNodeMap();
175  std::vector<unsigned> indices2 = reader2.GetIncompleteNodeMap();
176 
177  if (indices1.size() != indices2.size())
178  {
179  std::cout << "Index map sizes " << indices1.size() << " and " << indices2.size() << " don't match\n";
180  return false;
181  }
182 
183  for (unsigned index=0; index<indices1.size(); index++)
184  {
185  if (indices1[index]!=indices2[index])
186  {
187  std::cout << "Time steps " << indices1[index] << " and " << indices2[index] << " don't match\n";
188  return false;
189  }
190  }
191 
192  // Check all the data
193  for (unsigned index=0; index<indices1.size(); index++)
194  {
195  unsigned node_index = indices1[index];
196  for (unsigned var=0; var<num_vars; var++)
197  {
198  std::vector<double> var_over_time1 = reader1.GetVariableOverTime(variable_names1[var], node_index);
199  std::vector<double> var_over_time2 = reader2.GetVariableOverTime(variable_names1[var], node_index);
200  for (unsigned time_step=0;time_step< var_over_time1.size(); time_step++)
201  {
202  if (fabs(var_over_time1[time_step] - var_over_time2[time_step]) > tol)
203  {
204  std::cout<<"Node "<<node_index<<" at time step "<<time_step<<" variable "<<variable_names1[var]<<
205  " differs ("<<var_over_time1[time_step]<<" != "<<var_over_time2[time_step]<<")\n";
206  return false;
207  }
208  }
209  }
210  }
211  }
212  return true;
213 }
214 
215 bool CompareFilesViaHdf5DataReaderGlobalNorm(std::string pathname1, std::string filename1, bool makeAbsolute1,
216  std::string pathname2, std::string filename2, bool makeAbsolute2,
217  double tol, std::string datasetName)
218 {
219  Hdf5DataReader reader1(pathname1, filename1, makeAbsolute1, datasetName);
220  Hdf5DataReader reader2(pathname2, filename2, makeAbsolute2, datasetName);
221 
222  unsigned number_nodes1 = reader1.GetNumberOfRows();
223  bool is_the_same = true;
224  std::vector<std::string> variable_names1 = reader1.GetVariableNames();
225  std::vector<std::string> variable_names2 = reader2.GetVariableNames();
226  std::vector<double> times1 = reader1.GetUnlimitedDimensionValues();
227  unsigned num_vars = variable_names1.size();
228  DistributedVectorFactory factory(number_nodes1);
229 
230  Vec data1 = factory.CreateVec();
231  Vec data2 = factory.CreateVec();
232 
233  for (unsigned timestep=0; timestep<times1.size(); timestep++)
234  {
235  for (unsigned var=0; var<num_vars; var++)
236  {
237  reader1.GetVariableOverNodes(data1, variable_names1[var], timestep);
238  reader2.GetVariableOverNodes(data2, variable_names2[var], timestep);
239 
240  PetscReal data1_norm;
241  PetscReal data2_norm;
242  VecNorm(data1, NORM_2, &data1_norm);
243  VecNorm(data2, NORM_2, &data2_norm);
244  PetscReal norm = fabs(data1_norm-data2_norm);
245  if (norm > tol)
246  {
247  is_the_same = false;
248  std::cout << "Vectors differ in global NORM_2 by " << norm << std::endl;
249  }
250  }
251  }
252 
253  PetscTools::Destroy(data1);
254  PetscTools::Destroy(data2);
255 
256  return is_the_same;
257 }
static void Destroy(Vec &rVec)
Definition: PetscTools.hpp:351