| 1 | |
|---|
| 2 | #ifndef __FIELD_IMPL_HPP__ |
|---|
| 3 | #define __FIELD_IMPL_HPP__ |
|---|
| 4 | |
|---|
| 5 | #include "xios_spl.hpp" |
|---|
| 6 | #include "field.hpp" |
|---|
| 7 | #include "context.hpp" |
|---|
| 8 | #include "grid.hpp" |
|---|
| 9 | #include "timer.hpp" |
|---|
| 10 | #include "array_new.hpp" |
|---|
| 11 | #include "source_filter.hpp" |
|---|
| 12 | #include "store_filter.hpp" |
|---|
| 13 | |
|---|
| 14 | |
|---|
| 15 | namespace xios { |
|---|
| 16 | |
|---|
| 17 | template <int N> |
|---|
| 18 | void CField::setData(const CArray<double, N>& _data, int tileid) |
|---|
| 19 | TRY |
|---|
| 20 | { |
|---|
| 21 | if (clientSourceFilter) |
|---|
| 22 | { |
|---|
| 23 | if (check_if_active.isEmpty() || (!check_if_active.isEmpty() && (!check_if_active) || isActive(true))) |
|---|
| 24 | { |
|---|
| 25 | if ( info.getLevel()>100 ) |
|---|
| 26 | { |
|---|
| 27 | const double* array = _data.dataFirst(); |
|---|
| 28 | int numElements( _data.numElements() ); |
|---|
| 29 | checkSumLike( array, numElements, true ); |
|---|
| 30 | } |
|---|
| 31 | if (tileid > -1) |
|---|
| 32 | clientSourceFilter->streamTile(CContext::getCurrent()->getCalendar()->getCurrentDate(), _data, tileid); // tiled domain |
|---|
| 33 | else |
|---|
| 34 | clientSourceFilter->streamData(CContext::getCurrent()->getCalendar()->getCurrentDate(), _data); |
|---|
| 35 | } |
|---|
| 36 | } |
|---|
| 37 | else if (instantDataFilter) |
|---|
| 38 | ERROR("void CField::setData(const CArray<double, N>& _data)", |
|---|
| 39 | << "Impossible to receive data from the model for a field [ id = " << getId() << " ] with a reference or an arithmetic operation."); |
|---|
| 40 | } |
|---|
| 41 | CATCH_DUMP_ATTR |
|---|
| 42 | |
|---|
| 43 | template <int N> |
|---|
| 44 | void CField::getData(CArray<double, N>& _data) const |
|---|
| 45 | TRY |
|---|
| 46 | { |
|---|
| 47 | if (storeFilter) |
|---|
| 48 | { |
|---|
| 49 | CDataPacket::StatusCode status = storeFilter->getData(CContext::getCurrent()->getCalendar()->getCurrentDate(), _data); |
|---|
| 50 | if ( info.getLevel()>100 ) |
|---|
| 51 | { |
|---|
| 52 | const double* array = _data.dataFirst(); |
|---|
| 53 | int numElements( _data.numElements() ); |
|---|
| 54 | checkSumLike( array, numElements, false ); |
|---|
| 55 | } |
|---|
| 56 | |
|---|
| 57 | if (status == CDataPacket::END_OF_STREAM) |
|---|
| 58 | ERROR("void CField::getData(CArray<double, N>& _data) const", |
|---|
| 59 | << "Impossible to access field data, all the records of the field [ id = " << getId() << " ] have been already read."); |
|---|
| 60 | } |
|---|
| 61 | else |
|---|
| 62 | { |
|---|
| 63 | ERROR("void CField::getData(CArray<double, N>& _data) const", |
|---|
| 64 | << "Impossible to access field data, the field [ id = " << getId() << " ] does not have read access."); |
|---|
| 65 | } |
|---|
| 66 | } |
|---|
| 67 | CATCH |
|---|
| 68 | |
|---|
| 69 | void CField::checkSumLike( const double* array, int numElements, bool output ) const |
|---|
| 70 | { |
|---|
| 71 | int rk = CContext::getCurrent()->client->clientRank; |
|---|
| 72 | int sz = CContext::getCurrent()->client->clientSize; |
|---|
| 73 | MPI_Comm comm = CContext::getCurrent()->client->intraComm; |
|---|
| 74 | |
|---|
| 75 | double localSum( 0. ); |
|---|
| 76 | double error( 0. ); |
|---|
| 77 | unsigned long long checkSum( 0 ); |
|---|
| 78 | for ( int i=0 ; i<numElements ; i++ ) { |
|---|
| 79 | bool contributes( true ); |
|---|
| 80 | if ( (!output) && ( !default_value.isEmpty() ) ) |
|---|
| 81 | { |
|---|
| 82 | if ( fabs(array[i]) > 0) |
|---|
| 83 | if ( fabs(array[i]-default_value.getValue())/array[i] < 2e-16 ) |
|---|
| 84 | contributes = false; |
|---|
| 85 | } |
|---|
| 86 | if (contributes) |
|---|
| 87 | { |
|---|
| 88 | double y = array[i] - error; |
|---|
| 89 | double t = localSum + y; |
|---|
| 90 | error = (t - localSum) - y; |
|---|
| 91 | localSum = t; |
|---|
| 92 | |
|---|
| 93 | checkSum += *((unsigned long long*)(&array[i])); |
|---|
| 94 | checkSum = checkSum%LLONG_MAX; |
|---|
| 95 | } |
|---|
| 96 | } |
|---|
| 97 | |
|---|
| 98 | double globalSum( 0. ); |
|---|
| 99 | unsigned long long globalCheck( 0 ); |
|---|
| 100 | |
|---|
| 101 | if ( rk ==0 ) |
|---|
| 102 | { |
|---|
| 103 | MPI_Status status; |
|---|
| 104 | globalSum = localSum; // rank 0 contribution |
|---|
| 105 | globalCheck = checkSum; |
|---|
| 106 | for ( int irk = 1 ; irk < sz ; irk++ ) |
|---|
| 107 | { |
|---|
| 108 | MPI_Recv( &localSum, 1, MPI_DOUBLE, irk, 0, comm, &status ); |
|---|
| 109 | double y = localSum - error; |
|---|
| 110 | double t = globalSum + y; |
|---|
| 111 | error = (t - globalSum) - y; |
|---|
| 112 | globalSum = t; |
|---|
| 113 | |
|---|
| 114 | MPI_Recv( &checkSum, 1, MPI_UNSIGNED_LONG_LONG, irk, 1, comm, &status ); |
|---|
| 115 | globalCheck += checkSum; |
|---|
| 116 | globalCheck = globalCheck%LLONG_MAX; |
|---|
| 117 | } |
|---|
| 118 | } |
|---|
| 119 | else |
|---|
| 120 | { |
|---|
| 121 | MPI_Send( &localSum, 1, MPI_DOUBLE, 0, 0, comm ); |
|---|
| 122 | MPI_Send( &checkSum, 1, MPI_UNSIGNED_LONG_LONG, 0, 1, comm ); |
|---|
| 123 | } |
|---|
| 124 | |
|---|
| 125 | if ( rk == 0 ) |
|---|
| 126 | { |
|---|
| 127 | info(100) << setprecision(DBL_DIG); |
|---|
| 128 | if (output ) |
|---|
| 129 | info(100) << "Check Output key field for : " << getId() << ", key = " << globalCheck << ", sum = " << globalSum << endl; |
|---|
| 130 | else |
|---|
| 131 | info(100) << "Check Input key field for : " << getId() << ", key = " << globalCheck << ", sum = " << globalSum << endl; |
|---|
| 132 | info(100) << setprecision(6); |
|---|
| 133 | } |
|---|
| 134 | } |
|---|
| 135 | |
|---|
| 136 | } // namespace xios |
|---|
| 137 | |
|---|
| 138 | #endif |
|---|
