/ return L2 norm of error normalized by average of L2 norms */ /
| 103 | /* return L2 norm of error normalized by average of L2 norms */ |
| 104 | /***************************************************************/ |
| 105 | double compare_array_to_dataset(std::complex<meep::realnum> *field_array, int array_rank, |
| 106 | size_t *array_dims, const char *file, const char *name) { |
| 107 | int file_rank; |
| 108 | size_t file_dims[3]; |
| 109 | h5file f(file, h5file::READONLY, false); |
| 110 | char dataname[100]; |
| 111 | snprintf(dataname, 100, "%s.r", name); |
| 112 | double *rdata = |
| 113 | (double *)f.read(dataname, &file_rank, file_dims, 2, false /* single_precision */); |
| 114 | snprintf(dataname, 100, "%s.i", name); |
| 115 | double *idata = |
| 116 | (double *)f.read(dataname, &file_rank, file_dims, 2, false /* single_precision */); |
| 117 | if (!rdata || !idata) return -1.0; |
| 118 | if (file_rank != array_rank) return -1.0; |
| 119 | for (int n = 0; n < file_rank; n++) |
| 120 | if (file_dims[n] != array_dims[n]) return -1.0; |
| 121 | |
| 122 | double NormArray = 0.0, NormFile = 0.0, NormDelta = 0.0; |
| 123 | for (size_t n = 0; n < file_dims[0] * file_dims[1]; n++) { |
| 124 | std::complex<double> zArray = field_array[n]; |
| 125 | std::complex<double> zFile = std::complex<double>(rdata[n], idata[n]); |
| 126 | NormArray += norm(zArray); |
| 127 | NormFile += norm(zFile); |
| 128 | NormDelta += norm(zArray - zFile); |
| 129 | } |
| 130 | NormArray = sqrt(NormArray); |
| 131 | NormFile = sqrt(NormFile); |
| 132 | NormDelta = sqrt(NormDelta); |
| 133 | double RelErr = NormDelta / (0.5 * (NormArray + NormFile)); |
| 134 | return RelErr; |
| 135 | } |
| 136 | |
| 137 | /***************************************************************/ |
| 138 | /* compute the L2 norm of two complex-valued HDF5 datasets, */ |