00001
00009 #ifndef BS_FI_OPERATOR_NORM_CALC_H_
00010 #define BS_FI_OPERATOR_NORM_CALC_H_
00011
00012 namespace blue_sky {
00013
00014 #define MAX_AND_INDEX(NEW,OLD,I,I_OLD) if (fabs (NEW) > fabs (OLD)) {(OLD) = (NEW); (I_OLD) = (I);}
00015 #ifdef _MPI
00016 struct di
00017 {
00018 double d;
00019 int i;
00020 };
00021 #define MPI_MAX_AND_INDEX(NEW,I,DI_OLD,DI_NEW) DI_OLD.d = fabs (NEW); DI_OLD.i = (NEW > 0) ? I : -I;\
00022 MPI_Allreduce (&DI_OLD, &DI_NEW, 1, MPI_DOUBLE_INT, MPI_MAXLOC, MPI_COMM_WORLD); \
00023 NEW = (DI_NEW.i > 0) ? DI_NEW.d : - DI_NEW.d; I = abs(DI_NEW.i);
00024 #endif //_MPI
00025
00029 template <typename strategy_t, bool is_w, bool is_g, bool is_o>
00030 BS_FORCE_INLINE void
00031 fi_operator_impl <strategy_t, is_w, is_g, is_o>::norm_calc ()
00032 {
00033 double pv = 0;
00034
00035 #ifdef _MPI
00036 int n_left = 0;
00037 int n_right = n_cells_;
00038 int n_offset = 0;
00039 double mpi_l2_norm[6], mpi_mat_balance[3];
00040 struct di max_ind_old, max_ind_new;
00041 #else //_MPI
00042 int n_left = 0;
00043 int n_right = n_cells_;
00044 #endif //_MPI
00045
00046 norm_.clear ();
00047 for (index_t i = n_left; i < n_right; ++i)
00048 {
00049
00050 pv += volume_[i] * data_[i].porosity;
00051 update_norm_by_cell (i, norm_);
00052 }
00053
00054 #ifdef _MPI
00055 BS_ASSERT (false && "MPI: NOT IMPL YET");
00056
00057 double mpi_pv;
00058 MPI_Allreduce (&pv, &mpi_pv, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
00059 pv = mpi_pv;
00060
00062
00063
00064 MPI_Allreduce (&norm.val[norms::L2_CPV_WATER], &mpi_l2_norm[0], 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
00065 MPI_Allreduce (&norm.val[norms::L2_CPV_GAS], &mpi_l2_norm[1], 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
00066 MPI_Allreduce (&norm.val[norms::L2_CPV_OIL], &mpi_l2_norm[2], 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
00067 MPI_Allreduce (&norm.val[norms::L2_ACPV_WATER], &mpi_l2_norm[3], 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
00068 MPI_Allreduce (&norm.val[norms::L2_ACPV_GAS], &mpi_l2_norm[4], 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
00069 MPI_Allreduce (&norm.val[norms::L2_ACPV_OIL], &mpi_l2_norm[5], 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
00070 MPI_Allreduce (&norm.val[norms::MB_ERR_WATER], &mpi_mat_balance[0], 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
00071 MPI_Allreduce (&norm.val[norms::MB_ERR_GAS], &mpi_mat_balance[1], 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
00072 MPI_Allreduce (&norm.val[norms::MB_ERR_OIL], &mpi_mat_balance[2], 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
00073
00074
00075 norm.val[norms::L2_CPV_WATER] = mpi_l2_norm[0];
00076 norm.val[norms::L2_CPV_GAS] = mpi_l2_norm[1];
00077 norm.val[norms::L2_CPV_OIL] = mpi_l2_norm[2];
00078 norm.val[norms::L2_ACPV_WATER] = mpi_l2_norm[3];
00079 norm.val[norms::L2_ACPV_GAS] = mpi_l2_norm[4];
00080 norm.val[norms::L2_ACPV_OIL] = mpi_l2_norm[5];
00081 norm.val[norms::MB_ERR_WATER] = mpi_mat_balance[0];
00082 norm.val[norms::MB_ERR_GAS] = mpi_mat_balance[1];
00083 norm.val[norms::MB_ERR_OIL] = mpi_mat_balance[2];
00084
00085 MPI_MAX_AND_INDEX (norm.val[norms::C_CPV_WATER], norm.idx[norms::C_CPV_WATER], max_ind_old, max_ind_new);
00086 MPI_MAX_AND_INDEX (norm.val[norms::C_CPV_GAS], norm.idx[norms::C_CPV_GAS], max_ind_old, max_ind_new);
00087 MPI_MAX_AND_INDEX (norm.val[norms::C_CPV_OIL], norm.idx[norms::C_CPV_OIL], max_ind_old, max_ind_new);
00088 MPI_MAX_AND_INDEX (norm.val[norms::C_ACPV_WATER], norm.idx[norms::C_ACPV_WATER], max_ind_old, max_ind_new);
00089 MPI_MAX_AND_INDEX (norm.val[norms::C_ACPV_GAS], norm.idx[norms::C_ACPV_GAS], max_ind_old, max_ind_new);
00090 MPI_MAX_AND_INDEX (norm.val[norms::C_ACPV_OIL], norm.idx[norms::C_ACPV_OIL], max_ind_old, max_ind_new);
00091
00093 #endif //_MPI
00094
00095 double n = 1.0 / double (n_cells_);
00096 norm_.val[norms::L2_CPV_WATER] = sqrt (norm_.val[norms::L2_CPV_WATER]) * n;
00097 norm_.val[norms::L2_CPV_GAS] = sqrt (norm_.val[norms::L2_CPV_GAS]) * n;
00098 norm_.val[norms::L2_CPV_OIL] = sqrt (norm_.val[norms::L2_CPV_OIL]) * n;
00099 norm_.val[norms::L2_ACPV_WATER] = sqrt (norm_.val[norms::L2_ACPV_WATER]) * n;
00100 norm_.val[norms::L2_ACPV_GAS] = sqrt (norm_.val[norms::L2_ACPV_GAS]) * n;
00101 norm_.val[norms::L2_ACPV_OIL] = sqrt (norm_.val[norms::L2_ACPV_OIL]) * n;
00102
00103
00104 if (is_w)
00105 {
00106 MAX_AND_INDEX (norm_.val[norms::C_CPV_WATER], norm_.val[norms::C_CPV],
00107 norm_.idx[norms::C_CPV_WATER], norm_.idx[norms::C_CPV]);
00108 norm_.val[norms::C_ACPV_WATER] /= calc_model_->ave_volume;
00109
00110 MAX_AND_INDEX (norm_.val[norms::C_ACPV_WATER], norm_.val[norms::C_ACPV],
00111 norm_.idx[norms::C_ACPV_WATER], norm_.idx[norms::C_ACPV]);
00112
00113 MAX_AND_INDEX (norm_.val[norms::L2_CPV_WATER], norm_.val[norms::L2_CPV],
00114 norm_.idx[norms::L2_CPV_WATER], norm_.idx[norms::L2_CPV]);
00115
00116 MAX_AND_INDEX (norm_.val[norms::L2_ACPV_WATER], norm_.val[norms::L2_ACPV],
00117 norm_.idx[norms::L2_ACPV_WATER], norm_.idx[norms::L2_ACPV]);
00118
00119 norm_.val[norms::MB_ERR_WATER] /= pv * calc_model_->invers_fvf_average[d_w];
00120 MAX_AND_INDEX (norm_.val[norms::MB_ERR_WATER], norm_.val[norms::MB_ERR],
00121 norm_.idx[norms::MB_ERR_WATER], norm_.idx[norms::MB_ERR]);
00122 }
00123 if (is_o)
00124 {
00125 MAX_AND_INDEX (norm_.val[norms::C_CPV_OIL], norm_.val[norms::C_CPV],
00126 norm_.idx[norms::C_CPV_OIL], norm_.idx[norms::C_CPV]);
00127 norm_.val[norms::C_ACPV_OIL] /= calc_model_->ave_volume;
00128
00129 MAX_AND_INDEX (norm_.val[norms::C_ACPV_OIL], norm_.val[norms::C_ACPV],
00130 norm_.idx[norms::C_ACPV_OIL], norm_.idx[norms::C_ACPV]);
00131
00132 MAX_AND_INDEX (norm_.val[norms::L2_CPV_OIL], norm_.val[norms::L2_CPV],
00133 norm_.idx[norms::L2_CPV_OIL], norm_.idx[norms::L2_CPV]);
00134
00135 MAX_AND_INDEX (norm_.val[norms::L2_ACPV_OIL], norm_.val[norms::L2_ACPV],
00136 norm_.idx[norms::L2_ACPV_OIL], norm_.idx[norms::L2_ACPV]);
00137
00138 norm_.val[norms::MB_ERR_OIL] /= pv * calc_model_->invers_fvf_average[d_o];
00139 MAX_AND_INDEX (norm_.val[norms::MB_ERR_OIL], norm_.val[norms::MB_ERR],
00140 norm_.idx[norms::MB_ERR_OIL], norm_.idx[norms::MB_ERR]);
00141
00142 }
00143 if (is_g)
00144 {
00145 MAX_AND_INDEX (norm_.val[norms::C_CPV_GAS], norm_.val[norms::C_CPV],
00146 norm_.idx[norms::C_CPV_GAS], norm_.idx[norms::C_CPV]);
00147 norm_.val[norms::C_ACPV_GAS] /= calc_model_->ave_volume;
00148
00149 MAX_AND_INDEX (norm_.val[norms::C_ACPV_GAS], norm_.val[norms::C_ACPV],
00150 norm_.idx[norms::C_ACPV_GAS], norm_.idx[norms::C_ACPV]);
00151
00152 MAX_AND_INDEX (norm_.val[norms::L2_CPV_GAS], norm_.val[norms::L2_CPV],
00153 norm_.idx[norms::L2_CPV_GAS], norm_.idx[norms::L2_CPV]);
00154
00155 MAX_AND_INDEX (norm_.val[norms::L2_ACPV_GAS], norm_.val[norms::L2_ACPV],
00156 norm_.idx[norms::L2_ACPV_GAS], norm_.idx[norms::L2_ACPV]);
00157
00158 norm_.val[norms::MB_ERR_GAS] /= pv * calc_model_->invers_fvf_average[d_g];
00159 MAX_AND_INDEX (norm_.val[norms::MB_ERR_GAS], norm_.val[norms::MB_ERR],
00160 norm_.idx[norms::MB_ERR_GAS], norm_.idx[norms::MB_ERR]);
00161 }
00162
00163
00164 #ifdef _MPI
00165 if ((norm_.idx[norms::C_CPV] - n_offset >= n_left) && (norm_.idx[norms::C_CPV] - n_offset < n_right))
00166 {
00167 calc_model_->max_norm_counter[norm_.idx[norms::C_CPV] - n_offset]++;
00168 }
00169 #else //_MPI
00170 calc_model_->max_norm_counter[norm_.idx[norms::C_CPV]]++;
00171 #endif //_MPI
00172 }
00173
00179 template <class strategy_t, bool is_w, bool is_g, bool is_o>
00180 BS_FORCE_INLINE void
00181 fi_operator_impl <strategy_t, is_w, is_g, is_o>::update_norm_by_cell (index_t i, norms_storage_t &ns)
00182 {
00183 double norm_mult, r, d;
00184
00185 int equ_w, equ_g, equ_o;
00186 if (n_phases == 3)
00187 {
00188 equ_w = i * n_phases + p3_wat;
00189 equ_g = i * n_phases + p3_gas;
00190 equ_o = i * n_phases + p3_oil;
00191 }
00192 else if (n_phases == 2 && is_w)
00193 {
00194 equ_w = i * n_phases + p2ow_wat;
00195 equ_o = i * n_phases + p2ow_oil;
00196 }
00197 else if (n_phases == 2 && is_g)
00198 {
00199 equ_g = i * n_phases + p2og_gas;
00200 equ_o = i * n_phases + p2og_oil;
00201 }
00202 else
00203 {
00204 equ_g = equ_w = equ_o = i * n_phases + 0;
00205 }
00206 #ifdef _MPI
00207 int n_left = 0;
00208 int n_offset = 0;
00209
00210
00211 #endif //_MPI
00212
00213
00214 norm_mult = 1.0 / (volume_[i] * data_[i].porosity);
00215
00216 #ifdef _MPI
00217
00218 i += n_offset;
00219 #endif //_MPI
00220
00221
00222 if (is_w)
00223 {
00224 r = rhs_[equ_w] + flux_rhs_[equ_w];
00225
00226
00227 ns.val[norms::MB_ERR_WATER] += r / calc_model_->invers_fvf_average[d_w];
00228
00229
00230 d = r * norm_mult / calc_model_->invers_fvf_average[d_w];
00231 MAX_AND_INDEX (d, ns.val[norms::C_CPV_WATER], i, ns.idx[norms::C_CPV_WATER]);
00232
00233
00234 ns.val[norms::L2_CPV_WATER] += d * d;
00235
00236
00237 d = r / data_[i].invers_fvf[d_w];
00238 MAX_AND_INDEX (d, ns.val[norms::C_ACPV_WATER], i, ns.idx[norms::C_ACPV_WATER]);
00239
00240
00241 ns.val[norms::L2_ACPV_WATER] += d * d;
00242 }
00243
00244 if (is_g)
00245 {
00246 r = rhs_[equ_g] + flux_rhs_[equ_g];
00247
00248
00249 ns.val[norms::MB_ERR_GAS] += r / calc_model_->invers_fvf_average[d_g];
00250
00251
00252 d = r * norm_mult / calc_model_->invers_fvf_average[d_g];
00253 MAX_AND_INDEX (d, ns.val[norms::C_CPV_GAS], i, ns.idx[norms::C_CPV_GAS]);
00254
00255
00256 ns.val[norms::L2_CPV_GAS] += d * d;
00257
00258
00259 d = r / data_[i].invers_fvf[d_g];
00260 MAX_AND_INDEX (d, ns.val[norms::C_ACPV_GAS], i, ns.idx[norms::C_ACPV_GAS]);
00261
00262
00263 ns.val[norms::L2_ACPV_GAS] += d * d;
00264 }
00265
00266 if (is_o)
00267 {
00268 r = rhs_[equ_o] + flux_rhs_[equ_o];
00269
00270
00271 ns.val[norms::MB_ERR_OIL] += r / calc_model_->invers_fvf_average[d_o];
00272
00273
00274 d = r * norm_mult / calc_model_->invers_fvf_average[d_o];
00275 MAX_AND_INDEX (d, ns.val[norms::C_CPV_OIL], i, ns.idx[norms::C_CPV_OIL]);
00276
00277
00278 ns.val[norms::L2_CPV_OIL] += d * d;
00279
00280
00281 d = r / data_[i].invers_fvf[d_o];
00282 MAX_AND_INDEX (d, ns.val[norms::C_ACPV_OIL], i, ns.idx[norms::C_ACPV_OIL]);
00283
00284
00285 ns.val[norms::L2_ACPV_OIL] += d * d;
00286 }
00287
00288 rhs_item_array_t &s_rhs = jmatrix_->get_sec_rhs ();
00289 for (index_t j = 0, j_cnt = n_sec_vars; j < j_cnt; ++j)
00290 {
00291 MAX_AND_INDEX (s_rhs[i * n_sec_vars + j], ns.val[norms::S_RHS], i, ns.idx[norms::S_RHS]);
00292 }
00293 }
00294
00295 }
00296
00297 #endif // #ifndef BS_FI_OPERATOR_NORM_CALC_H_