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 | #include <cfloat> |
---|
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 ( CXios::checkSumSend ) |
---|
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 ( CXios::checkSumRecv ) |
---|
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 |
---|