CompareHdf5ResultsFiles.cpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2011
00004 
00005 University of Oxford means the Chancellor, Masters and Scholars of the
00006 University of Oxford, having an administrative office at Wellington
00007 Square, Oxford OX1 2JD, UK.
00008 
00009 This file is part of Chaste.
00010 
00011 Chaste is free software: you can redistribute it and/or modify it
00012 under the terms of the GNU Lesser General Public License as published
00013 by the Free Software Foundation, either version 2.1 of the License, or
00014 (at your option) any later version.
00015 
00016 Chaste is distributed in the hope that it will be useful, but WITHOUT
00017 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
00019 License for more details. The offer of Chaste under the terms of the
00020 License is subject to the License being interpreted in accordance with
00021 English Law and subject to any action against the University of Oxford
00022 being under the jurisdiction of the English Courts.
00023 
00024 You should have received a copy of the GNU Lesser General Public License
00025 along with Chaste. If not, see <http://www.gnu.org/licenses/>.
00026 
00027 */
00028 
00029 #include "CompareHdf5ResultsFiles.hpp"
00030 
00031 #include <iostream>
00032 #include <vector>
00033 #include "petscvec.h"
00034 #include "Hdf5DataReader.hpp"
00035 #include "DistributedVectorFactory.hpp"
00036 
00037 
00038 bool CompareFilesViaHdf5DataReader(std::string pathname1, std::string filename1, bool makeAbsolute1,
00039                                    std::string pathname2, std::string filename2, bool makeAbsolute2,
00040                                    double tol)
00041 {
00042     Hdf5DataReader reader1(pathname1, filename1, makeAbsolute1);
00043     Hdf5DataReader reader2(pathname2, filename2, makeAbsolute2);
00044 
00045     unsigned number_nodes1 = reader1.GetNumberOfRows();
00046     unsigned number_nodes2 = reader2.GetNumberOfRows();
00047     if (number_nodes1 != number_nodes2)
00048     {
00049         std::cout << "Number of nodes " << number_nodes1 << " and " << number_nodes2 << " don't match\n";
00050         return false;
00051     }
00052     // Check the variable names and units
00053     std::vector<std::string> variable_names1 = reader1.GetVariableNames();
00054     std::vector<std::string> variable_names2 = reader2.GetVariableNames();
00055     unsigned num_vars = variable_names1.size();
00056     if (num_vars != variable_names2.size())
00057     {
00058         std::cout << "Number of variables " << variable_names1.size()
00059                   << " and " << variable_names2.size() << " don't match\n";
00060         return false;
00061     }
00062     for (unsigned var=0; var<num_vars; var++)
00063     {
00064         std::string var_name = variable_names1[var];
00065         if (var_name != variable_names2[var])
00066         {
00067             std::cout << "Variable names " << var_name << " and "
00068                       << variable_names2[var] << " don't match\n";
00069             return false;
00070         }
00071         if (reader1.GetUnit(var_name) != reader2.GetUnit(var_name))
00072         {
00073             std::cout << "Units names " << reader1.GetUnit(var_name)
00074                       << " and " << reader2.GetUnit(var_name) << " don't match\n";
00075             return false;
00076         }
00077     }
00078     // Check the timestep vectors
00079     std::vector<double> times1 = reader1.GetUnlimitedDimensionValues();
00080     std::vector<double> times2 = reader2.GetUnlimitedDimensionValues();
00081 
00082     if (times1.size() != times2.size())
00083     {
00084         std::cout << "Time step sizes " << times1.size()
00085                   << " and " << times2.size() << " don't match\n";
00086         return false;
00087     }
00088 
00089     for (unsigned timestep=0; timestep<times1.size(); timestep++)
00090     {
00091         if (fabs(times1[timestep]-times2[timestep]) > 1e-8)
00092         {
00093             std::cout << "Time steps " << times1[timestep]
00094                       << " and " << times2[timestep] << " don't match\n";
00095             return false;
00096         }
00097     }
00098 
00099     bool is_complete1 = reader1.IsDataComplete();
00100     bool is_complete2 = reader2.IsDataComplete();
00101 
00102     if (is_complete1 != is_complete2)
00103     {
00104         std::cout<<"One of the readers has incomplete data and the other doesn't\n";
00105         return false;
00106     }
00107 
00108     if (is_complete1)
00109     {
00110         DistributedVectorFactory factory(number_nodes1);
00111 
00112         Vec data1 = factory.CreateVec();
00113         Vec data2 = factory.CreateVec();
00114 
00115         for (unsigned timestep=0; timestep<times1.size(); timestep++)
00116         {
00117             for (unsigned var=0; var<num_vars; var++)
00118             {
00119                 reader1.GetVariableOverNodes(data1, variable_names1[var], timestep);
00120                 reader2.GetVariableOverNodes(data2, variable_names2[var], timestep);
00121 
00122 #if (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 2) //PETSc 2.2
00123                 double minus_one = -1.0;
00124                 VecAXPY(&minus_one, data2, data1);
00125 #else
00126                 //[note: VecAXPY(y,a,x) computes y = ax+y]
00127                 VecAXPY(data1, -1.0, data2);
00128 #endif
00129 
00130                 PetscReal difference_norm;
00131                 VecNorm(data1, NORM_2, &difference_norm);
00132 
00133                 if (difference_norm > tol)
00134                 {
00135                     std::cout << "Vectors differ in NORM_2 by " << difference_norm << std::endl;
00136                     return false;
00137                 }
00138             }
00139         }
00140         VecDestroy(data1);
00141         VecDestroy(data2);
00142     }
00143     else
00144     {
00145         // Incomplete data
00146 
00147         // Check the index vectors
00148         std::vector<unsigned> indices1 = reader1.GetIncompleteNodeMap();
00149         std::vector<unsigned> indices2 = reader2.GetIncompleteNodeMap();
00150 
00151         if (indices1.size() != indices2.size())
00152         {
00153             std::cout << "Index map sizes " << indices1.size() << " and " << indices2.size() << " don't match\n";
00154             return false;
00155         }
00156 
00157         for (unsigned index=0; index<indices1.size(); index++)
00158         {
00159             if (indices1[index]!=indices2[index])
00160             {
00161                std::cout << "Time steps " << indices1[index] << " and " << indices2[index] << " don't match\n";
00162                return false;
00163             }
00164         }
00165 
00166         // Check all the data
00167         for (unsigned index=0; index<indices1.size(); index++)
00168         {
00169             unsigned node_index = indices1[index];
00170             for (unsigned var=0; var<num_vars; var++)
00171             {
00172                 std::vector<double> var_over_time1 = reader1.GetVariableOverTime(variable_names1[var], node_index);
00173                 std::vector<double> var_over_time2 = reader2.GetVariableOverTime(variable_names1[var], node_index);
00174                 for (unsigned time_step=0;time_step< var_over_time1.size(); time_step++)
00175                 {
00176                     if (fabs(var_over_time1[time_step] - var_over_time2[time_step]) > tol)
00177                     {
00178                         std::cout<<"Node "<<node_index<<" at time step "<<time_step<<" variable "<<variable_names1[var]<<
00179                             " differs ("<<var_over_time1[time_step]<<" != "<<var_over_time2[time_step]<<")\n";
00180                         return false;
00181                     }
00182                 }
00183             }
00184         }
00185     }
00186     return true;
00187 }
00188 
00189 bool CompareFilesViaHdf5DataReaderGlobalNorm(std::string pathname1, std::string filename1, bool makeAbsolute1,
00190                                              std::string pathname2, std::string filename2, bool makeAbsolute2,
00191                                              double tol)
00192 {
00193         Hdf5DataReader reader1(pathname1, filename1, makeAbsolute1);
00194         Hdf5DataReader reader2(pathname2, filename2, makeAbsolute2);
00195 
00196         unsigned number_nodes1 = reader1.GetNumberOfRows();
00197         bool is_the_same = true;
00198         std::vector<std::string> variable_names1 = reader1.GetVariableNames();
00199         std::vector<std::string> variable_names2 = reader2.GetVariableNames();
00200         std::vector<double> times1 = reader1.GetUnlimitedDimensionValues();
00201         unsigned num_vars = variable_names1.size();
00202         DistributedVectorFactory factory(number_nodes1);
00203 
00204         Vec data1 = factory.CreateVec();
00205         Vec data2 = factory.CreateVec();
00206 
00207         for (unsigned timestep=0; timestep<times1.size(); timestep++)
00208         {
00209             for (unsigned var=0; var<num_vars; var++)
00210             {
00211                 reader1.GetVariableOverNodes(data1, variable_names1[var], timestep);
00212                 reader2.GetVariableOverNodes(data2, variable_names2[var], timestep);
00213 
00214                 PetscReal data1_norm;
00215                 PetscReal data2_norm;
00216                 VecNorm(data1, NORM_2, &data1_norm);
00217                 VecNorm(data2, NORM_2, &data2_norm);
00218                 PetscReal norm = fabs(data1_norm-data2_norm);
00219                 if (norm > tol)
00220                 {
00221                     is_the_same = false;
00222                     std::cout << "Vectors differ in global NORM_2 by " << norm << std::endl;
00223                 }
00224             }
00225         }
00226 
00227         VecDestroy(data1);
00228         VecDestroy(data2);
00229 
00230         return is_the_same;
00231 }

Generated on Mon Apr 18 11:35:28 2011 for Chaste by  doxygen 1.5.5