00001
00010 #ifndef BS_FI_OPERATOR_H_
00011 #define BS_FI_OPERATOR_H_
00012
00013 #include BS_FORCE_PLUGIN_IMPORT ()
00014 #include "switch_main_vars.h"
00015 #include "pvt_water.h"
00016 #include "pvt_dead_oil.h"
00017 #include "pvt_gas.h"
00018 #include "scal_region_info.h"
00019 #include "scal_region.h"
00020 #include "scal_2p_data_holder.h"
00021 #include "scal_3p.h"
00022 #include BS_STOP_PLUGIN_IMPORT ()
00023
00024 #include "calc_model.h"
00025 #include "norm_calc.h"
00026
00027 namespace blue_sky {
00028
00035 template <typename strategy_type, bool is_w, bool is_g, bool is_o>
00036 struct fi_operator_impl
00037 {
00038 typedef strategy_type strategy_t;
00039 typedef typename strategy_t::item_t item_t;
00040 typedef typename strategy_t::rhs_item_t rhs_item_t;
00041 typedef typename strategy_t::index_t index_t;
00042 typedef typename strategy_t::item_array_t item_array_t;
00043 typedef typename strategy_t::rhs_item_array_t rhs_item_array_t;
00044 typedef typename strategy_t::index_array_t index_array_t;
00045 typedef typename strategy_t::csr_matrix_t bcsr_matrix_t;
00046
00047 typedef calc_model <strategy_t> calc_model_t;
00048 typedef typename calc_model_t::data_t data_t;
00049 typedef typename calc_model_t::data_array_t data_array_t;
00050 typedef typename calc_model_t::main_var_array_t main_var_array_t;
00051 typedef norms_storage <strategy_t> norms_storage_t;
00052 typedef reservoir <strategy_t> reservoir_t;
00053 typedef rs_mesh_iface <strategy_t> mesh_iface_t;
00054 typedef jacobian <strategy_t> jacobian_t;
00055 typedef jacobian_matrix <strategy_t> jmatrix_t;
00056
00057 typedef typename calc_model_t::sp_this_t sp_calc_model_t;
00058 typedef typename calc_model_t::sp_reservoir_t sp_reservoir_t;
00059 typedef typename calc_model_t::sp_jacobian_t sp_jacobian_t;
00060 typedef typename calc_model_t::sp_jacobian_matrix_t sp_jmatrix_t;
00061 typedef typename calc_model_t::sp_mesh_iface_t sp_mesh_iface_t;
00062 typedef typename calc_model_t::sp_rock_grid sp_rock_grid_prop_t;
00063 typedef smart_ptr <bcsr_matrix_t, true> sp_bcsr_matrix_t;
00064
00065 typedef typename calc_model_t::sp_pvt_dead_oil_array_t sp_pvt_dead_oil_array_t;
00066 typedef typename calc_model_t::sp_pvt_gas_array_t sp_pvt_gas_array_t;
00067 typedef typename calc_model_t::sp_pvt_water_array_t sp_pvt_water_array_t;
00068
00069 enum {
00070 n_phases = is_w + is_g + is_o,
00071 b_sqr = n_phases * n_phases,
00072 is_1p = n_phases == 1,
00073 };
00074
00075 public:
00084 fi_operator_impl (sp_calc_model_t &calc_model, sp_reservoir_t &reservoir, const sp_mesh_iface_t &mesh, sp_jacobian_t &jacobian, sp_jmatrix_t &jmatrix)
00085 : calc_model_ (calc_model),
00086 rock_grid_prop_ (calc_model_->rock_grid_prop),
00087 reservoir_ (reservoir),
00088 jacobian_ (jacobian),
00089 jmatrix_ (jmatrix),
00090 mesh_ (mesh),
00091 trns_matrix_ (jmatrix_->trns_matrix),
00092 trns_values_ (trns_matrix_->get_values ()),
00093 trns_rows_ptr_ (trns_matrix_->get_rows_ptr ()),
00094 trns_cols_ptr_ (trns_matrix_->get_cols_ind ()),
00095 reg_matrix_ (jmatrix_->get_regular_matrix ()),
00096 reg_values_ (reg_matrix_->get_values ()),
00097 reg_rows_ptr_ (reg_matrix_->get_rows_ptr ()),
00098 reg_cols_ptr_ (reg_matrix_->get_cols_ind ()),
00099 m_array_ (jmatrix_->m_array),
00100 p_array_ (jmatrix_->p_array),
00101 rhs_ (jmatrix_->get_rhs ()),
00102 sol_ (jmatrix_->get_solution ()),
00103 flux_rhs_ (jmatrix_->get_rhs_flux ()),
00104 sp_diag_ (jmatrix_->get_sp_diagonal ()),
00105 s_rhs_ (jmatrix_->get_sec_rhs ()),
00106 depths_ (mesh_->get_depths ()),
00107 n_cells_ (mesh->get_n_active_elements()),
00108 n_connections_ (mesh_->get_n_connections ()),
00109 n_sec_vars (calc_model_->n_sec_vars),
00110 d_w (calc_model_->phase_d[FI_PHASE_WATER]),
00111 d_g (calc_model_->phase_d[FI_PHASE_GAS]),
00112 d_o (calc_model_->phase_d[FI_PHASE_OIL]),
00113 ds_w (calc_model_->sat_d[FI_PHASE_WATER]),
00114 ds_g (calc_model_->sat_d[FI_PHASE_GAS]),
00115 d_gg (d_g * n_phases + d_g),
00116 d_gw (d_g * n_phases + d_w),
00117 d_go (d_g * n_phases + d_o),
00118 d_wg (d_w * n_phases + d_g),
00119 d_ww (d_w * n_phases + d_w),
00120 d_wo (d_w * n_phases + d_o),
00121 d_og (d_o * n_phases + d_g),
00122 d_ow (d_o * n_phases + d_w),
00123 d_oo (d_o * n_phases + d_o),
00124 data_ (calc_model_->data),
00125 saturation_3p_ (calc_model_->saturation_3p),
00126 pressure_ (calc_model_->pressure),
00127 gas_oil_ratio_ (calc_model_->gas_oil_ratio),
00128 main_vars_ (calc_model_->main_variable),
00129 volume_ (rock_grid_prop_->volume),
00130 poro_array_ (rock_grid_prop_->porosity_p_ref),
00131 rock_grid_comp_const_ (rock_grid_prop_->comp_const),
00132 rock_grid_comp_ref_pressure_ (rock_grid_prop_->comp_ref_pressure),
00133 sat_regions_ (calc_model_->sat_regions),
00134 pvt_regions_ (calc_model_->pvt_regions),
00135 pvt_oil_array (calc_model_->pvt_oil_array),
00136 pvt_water_array (calc_model_->pvt_water_array),
00137 pvt_gas_array (calc_model_->pvt_gas_array),
00138 min_p_ ((item_t) calc_model_->ts_params->get_float (fi_params::PVT_PRESSURE_RANGE_MIN)),
00139 max_p_ ((item_t) calc_model_->ts_params->get_float (fi_params::PVT_PRESSURE_RANGE_MAX)),
00140 drsdt_ ((item_t) calc_model_->ts_params->get_float (fi_params::DRSDT)),
00141 rhs_residual_ ((item_t) calc_model_->ts_params->get_float (fi_params::NEWTON_RESIDUAL)),
00142 mb_error_ ((item_t) calc_model_->ts_params->get_float (fi_params::MASS_BALANS_ERROR)),
00143 s_rhs_norm_ ((item_t) calc_model_->ts_params->get_float (fi_params::S_RHS_NORM)),
00144 norm_ (calc_model_->norm),
00145 cfl_ (jmatrix_->get_cfl_vector ())
00146 {
00147 }
00148
00160 fi_operator_return_type
00161 fi_operator (double &dt, index_t istart, index_t istart_well_contr, index_t &, bool update_rhs_after_gauss_elimination, bool save_debug_files)
00162 {
00163 save_debug_files;
00164 item_t mult = 1.0;
00165 bool tuning = calc_model_->ts_params->get_bool (fi_params::NEWTON_TUNING);
00166
00167 static norms_storage_t base_norm;
00168
00169 if (istart)
00170 {
00171 reservoir_->init_jacobian (jmatrix_, n_cells_);
00172 }
00173
00174 index_t n_approx = 5;
00175 for (index_t i = 0; i < n_approx; ++i)
00176 {
00177 mult = 1.0;
00178
00179
00180 jmatrix_->init (n_cells_, n_phases, 3, 0, n_sec_vars);
00181
00182 if (is_o && is_g)
00183 {
00184 fi_operator_switch_main_vars (dt);
00185 }
00186
00187 if (n_phases > 1)
00188 {
00189
00190 calc_model_->scal_prop->process (saturation_3p_,
00191 sat_regions_,
00192 rock_grid_prop_->permeability,
00193 poro_array_,
00194 data_);
00195 }
00196
00197
00198
00199 fi_operator_cells ((istart == 1) ? 1 : 0, dt);
00200 fi_operator_fill ();
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00213 if (block_connections_mpfa (dt))
00214 throw bs_exception("fi_operator_block_connections_mpfa", "return -1");
00215
00216 if (update_rhs_after_gauss_elimination)
00217 {
00218
00219 reservoir_->calc_wells (((istart_well_contr == 1) ? 1 : 0), dt, calc_model_, mesh_, jmatrix_);
00220
00221
00222 reservoir_->fill_rhs_wells (dt, calc_model_, flux_rhs_, update_rhs_after_gauss_elimination);
00223 }
00224
00225 norm_calc ();
00226
00227 if (update_rhs_after_gauss_elimination)
00228 {
00229
00230 reservoir_->end_jacobian (dt, calc_model_, jacobian_);
00231 }
00232
00233 #ifdef _DEBUG
00234
00235 {
00236 debug_save_data (dt);
00237 }
00238 #endif
00239
00240 if (check_norm (istart))
00241 return FI_OPERATOR_RETURN_FAIL;
00242
00243 if (check_solution_mult_cell (istart, i, n_approx, base_norm))
00244 continue;
00245
00246 do_dt_reduce (dt, istart);
00247 do_dt_tuning (dt, tuning);
00248
00249 break;
00250 }
00251
00252 base_norm = norm_;
00253
00254 jmatrix_->clear_solution ();
00255 jmatrix_->summ_rhs ();
00256
00257 return FI_OPERATOR_RETURN_OK;
00258 }
00259
00266 void
00267 fi_operator_init (index_t istart, double dt)
00268 {
00269 if (is_o && is_g)
00270 {
00271 fi_operator_switch_main_vars (dt);
00272 }
00273
00274 if (n_phases > 1)
00275 {
00276
00277 calc_model_->scal_prop->process (saturation_3p_,
00278 sat_regions_,
00279 rock_grid_prop_->permeability,
00280 poro_array_,
00281 data_);
00282 }
00283
00284 fi_operator_cells (istart, dt);
00285 }
00286
00292 void
00293 fi_operator_cells (index_t istart, const item_t dt)
00294 {
00295
00296 index_t i;
00297 rhs_item_t *jac_block = 0;
00298 rhs_item_t *rhs_block = 0;
00299 index_t b_sqr = n_phases * n_phases;
00300
00301 rhs_item_array_t &main_diag_acc = jmatrix_->get_regular_acc_diag();
00302 rhs_item_array_t &ss_diag = jmatrix_->get_ss_diagonal ();
00303 rhs_item_array_t &sp_diag = jmatrix_->get_sp_diagonal ();
00304 rhs_item_array_t &s_rhs = jmatrix_->get_sec_rhs ();
00305 rhs_item_array_t &rhs = jmatrix_->get_rhs ();
00306
00307 BS_ERROR (!main_diag_acc.empty (), "fi_operator_cells");
00308 BS_ERROR (!ss_diag.empty (), "fi_operator_cells");
00309 BS_ERROR (!sp_diag.empty (), "fi_operator_cells");
00310 BS_ERROR (!s_rhs.empty (), "fi_operator_cells");
00311 BS_ERROR (!rhs.empty (), "fi_operator_cells");
00312
00313 index_t switch_to_sg_count = 0;
00314 index_t switch_to_ro_count = 0;
00315 index_t switch_to_momg_count = 0;
00316 item_t total_volume = 0.0;
00317 rhs_item_t *ss_block = 0;
00318 rhs_item_t *sp_block = 0;
00319 rhs_item_t *s_rhs_block = 0;
00320
00321 #ifdef _MPI
00322 int n_left = 0;
00323 int n_own = 0;
00324
00325 double mpi_invers_fvf_average[3];
00326 int n_procs;
00327 #endif //_MPI
00328
00329 #ifdef FI_OPERATOR_CELLS_PARALLEL
00330 double invers_fvf_average_w = 0.;
00331 double invers_fvf_average_g = 0.;
00332 double invers_fvf_average_o = 0.;
00333 #endif //FI_OPERATOR_CELLS_PARALLEL
00334
00335
00336
00337 calc_model_->lsearch_force_newton_step = 0;
00338
00339
00340 assign (calc_model_->invers_fvf_average, 0);
00341
00342
00343 #ifdef FI_OPERATOR_CELLS_PARALLEL
00344 #pragma omp parallel for private (jac_block, rhs_block) \
00345 reduction (+:invers_fvf_average_w, invers_fvf_average_g, invers_fvf_average_o,total_volume)
00346 #endif //FI_OPERATOR_CELLS_PARALLEL
00347
00348 for (i = 0; i < n_cells_; ++i)
00349 {
00350 jac_block = &main_diag_acc[b_sqr * i];
00351 ss_block = &ss_diag[n_sec_vars * n_sec_vars * i];
00352 s_rhs_block = &s_rhs[n_sec_vars * i];
00353 rhs_block = &rhs[i * n_phases];
00354 #ifdef _MPI
00355 BS_ASSERT (false && "MPI: NOT IMPL YET");
00356 if ((i < n_left) || (i >= n_own))
00357 {
00360 }
00361 else
00362 {
00363 jac_block -= n_left * b_sqr;
00364 }
00365 #endif //_MPI
00366
00367
00368 #ifdef FI_OPERATOR_CELLS_PARALLEL
00369 fi_operator_cell (istart, dt, i, jac_block, rhs_block,
00370 ss_block, s_rhs_block,
00371 switch_to_sg_count, switch_to_ro_count, switch_to_momg_count,
00372 invers_fvf_average_w, invers_fvf_average_g, invers_fvf_average_o, total_volume);
00373 #else //FI_OPERATOR_CELLS_PARALLEL
00374 fi_operator_cell (istart, dt, i, jac_block, rhs_block,
00375 ss_block, s_rhs_block,
00376 switch_to_sg_count, switch_to_ro_count, switch_to_momg_count,
00377 total_volume);
00378 #endif //FI_OPERATOR_CELLS_PARALLEL
00379 }
00380
00381 #ifdef _MPI
00382 BS_ASSERT (false && "MPI: NOT IMPL YET");
00383 double mpi_total_volume;
00384 MPI_Allreduce (&total_volume, &mpi_total_volume, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
00385 total_volume = mpi_total_volume;
00386
00387 MPI_Allreduce (&invers_fvf_average, &mpi_invers_fvf_average, FI_PHASE_TOT, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
00388 if (is_w)
00389 invers_fvf_average[d_w] = mpi_invers_fvf_average[d_w];
00390 if (is_g)
00391 invers_fvf_average[d_g] = mpi_invers_fvf_average[d_g];
00392 if (is_o)
00393 invers_fvf_average[d_o] = mpi_invers_fvf_average[d_o];
00394
00396 #endif //_MPI
00397
00398 #ifdef FI_OPERATOR_CELLS_PARALLEL
00399 if (is_w)
00400 invers_fvf_average[d_w] = invers_fvf_average_w / (double) n_elements;
00401 if (is_g)
00402 invers_fvf_average[d_g] = invers_fvf_average_g / (double) n_elements;
00403 if (is_o)
00404 invers_fvf_average[d_o] = invers_fvf_average_o / (double) n_elements;
00405 #else //FI_OPERATOR_CELLS_PARALLEL
00406
00407 if (is_w)
00408 calc_model_->invers_fvf_average[d_w] /= (item_t) n_cells_;
00409 if (is_g)
00410 calc_model_->invers_fvf_average[d_g] /= (item_t) n_cells_;
00411 if (is_o)
00412 calc_model_->invers_fvf_average[d_o] /= (item_t) n_cells_;
00413 #endif //FI_OPERATOR_CELLS_PARALLEL
00414
00415 calc_model_->ave_volume = total_volume / (item_t) n_cells_;
00416 }
00417
00418 #define POROSITY data_i.porosity
00419 #define P_DERIV_POROSITY data_i.p_deriv_porosity
00420 #define INVERS_FVF_O data_i.invers_fvf[d_o]
00421 #define INVERS_FVF_G data_i.invers_fvf[d_g]
00422 #define INVERS_FVF_W data_i.invers_fvf[d_w]
00423 #define P_DERIV_INVERS_FVF_O data_i.p_deriv_invers_fvf[d_o]
00424 #define P_DERIV_INVERS_FVF_G data_i.p_deriv_invers_fvf[d_g]
00425 #define P_DERIV_INVERS_FVF_W data_i.p_deriv_invers_fvf[d_w]
00426 #define GAS_OIL_RATIO gas_oil_ratio_[i]
00427 #define P_DERIV_GAS_OIL_RATIO data_i.p_deriv_gas_oil_ratio
00428 #define VOLUME volume_[i]
00429 #define S_DERIV_CAP_PRESSURE_G data_i.s_deriv_cap_pressure[ds_g]
00430 #define S_DERIV_CAP_PRESSURE_W data_i.s_deriv_cap_pressure[ds_w]
00431 #define GOR_DERIV_INVERS_FVF data_i.gor_deriv_invers_fvf
00432 #define PREV_FLUID_VOLUME_O data_i.prev_fluid_volume[d_o]
00433 #define PREV_FLUID_VOLUME_G data_i.prev_fluid_volume[d_g]
00434 #define PREV_FLUID_VOLUME_W data_i.prev_fluid_volume[d_w]
00435
00440 struct local_data
00441 {
00442 item_t sat_w;
00443 item_t sat_g;
00444 item_t sat_o;
00445 boost::array <item_t, 3> ps_block;
00446 };
00447
00455 BS_FORCE_INLINE void
00456 fill_jacobian (index_t i, const data_t &data_i, rhs_item_t *jac_block, local_data &local_data_)
00457 {
00458
00459 if (is_w)
00460 {
00461 item_t p_der = VOLUME * (P_DERIV_POROSITY * INVERS_FVF_W + POROSITY * P_DERIV_INVERS_FVF_W);
00462 item_t sw_der = VOLUME * POROSITY * (INVERS_FVF_W + local_data_.sat_w * P_DERIV_INVERS_FVF_W * S_DERIV_CAP_PRESSURE_W);
00463
00464 if (n_phases == 3)
00465 {
00466 jac_block[p3_wat_po] = p_der * local_data_.sat_w;
00467 jac_block[p3_wat_sg] = 0;
00468 jac_block[p3_wat_so] = 0;
00469 local_data_.ps_block[p3_wat] = sw_der;
00470 }
00471 else if (n_phases == 2)
00472 {
00473 jac_block[p2ow_wat_po] = p_der * local_data_.sat_w;
00474 jac_block[p2ow_wat_so] = 0;
00475 local_data_.ps_block[p2ow_wat] = sw_der;
00476 }
00477 else
00478 {
00479 jac_block[0] = p_der;
00480 }
00481 }
00482
00483 if (is_g)
00484 {
00485 if (n_phases > 1)
00486 {
00487 item_t drg_p = 0, drg_sg = 0, drg_so = 0;
00488 if (FI_CHK_SG (main_vars_, i))
00489 {
00490
00491
00492
00493 drg_p = VOLUME * (local_data_.sat_g * (P_DERIV_POROSITY * INVERS_FVF_G + POROSITY * P_DERIV_INVERS_FVF_G) +
00494 local_data_.sat_o * (P_DERIV_POROSITY * INVERS_FVF_O * GAS_OIL_RATIO + POROSITY * (P_DERIV_GAS_OIL_RATIO * INVERS_FVF_O + GAS_OIL_RATIO * P_DERIV_INVERS_FVF_O)));
00495 drg_sg = VOLUME * POROSITY * (INVERS_FVF_G + local_data_.sat_g * P_DERIV_INVERS_FVF_G * S_DERIV_CAP_PRESSURE_G - GAS_OIL_RATIO * INVERS_FVF_O);
00496 }
00497 else if (FI_CHK_RO (main_vars_, i))
00498 {
00499 drg_p = VOLUME * local_data_.sat_o * GAS_OIL_RATIO * (P_DERIV_POROSITY * INVERS_FVF_O + P_DERIV_INVERS_FVF_O * POROSITY);
00500 drg_sg = VOLUME * POROSITY * local_data_.sat_o * (INVERS_FVF_O + GAS_OIL_RATIO * GOR_DERIV_INVERS_FVF);
00501 }
00502 drg_so = VOLUME * POROSITY * GAS_OIL_RATIO * INVERS_FVF_O;
00503 if (FI_CHK_MOMG (main_vars_, i))
00504 {
00505 drg_p = 0;
00506 drg_sg = VOLUME;
00507 drg_so = 0;
00508 }
00509 if (n_phases == 3)
00510 {
00511 jac_block[p3_gas_po] = drg_p;
00512 jac_block[p3_gas_sg] = drg_sg;
00513 jac_block[p3_gas_so] = drg_so;
00514 local_data_.ps_block[p3_gas] = 0;
00515 }
00516 else
00517 {
00518 jac_block[p2og_gas_po] = drg_p;
00519 jac_block[p2og_gas_sg] = drg_sg;
00520 local_data_.ps_block[p2og_gas] = drg_so;
00521 }
00522 }
00523 else
00524 {
00525 jac_block[0] = VOLUME * (P_DERIV_POROSITY * INVERS_FVF_G + P_DERIV_INVERS_FVF_G * POROSITY);
00526 }
00527
00528 }
00529
00530 if (is_o)
00531 {
00532 if (n_phases == 3)
00533 {
00534 jac_block[p3_oil_po] = VOLUME * local_data_.sat_o * (P_DERIV_POROSITY * INVERS_FVF_O + P_DERIV_INVERS_FVF_O * POROSITY);
00535 jac_block[p3_oil_so] = VOLUME * POROSITY * INVERS_FVF_O;
00536
00537 if (FI_CHK_SG (main_vars_, i))
00538 {
00539
00540 jac_block[p3_oil_sg] = -VOLUME * POROSITY * INVERS_FVF_O;
00541 }
00542 else if (FI_CHK_RO (main_vars_, i))
00543 {
00544 jac_block[p3_oil_sg] = VOLUME * POROSITY * local_data_.sat_o * GOR_DERIV_INVERS_FVF;
00545 }
00546 local_data_.ps_block[p3_oil] = 0.0;
00547
00548 if (FI_CHK_MOMG (main_vars_, i))
00549 {
00550 jac_block[p3_oil_po] = 0;
00551 jac_block[p3_oil_so] = VOLUME;
00552 jac_block[p3_oil_sg] = 0.0;
00553 }
00554 }
00555 else if (n_phases == 2 && is_g)
00556 {
00557 jac_block[p2og_oil_po] = VOLUME * local_data_.sat_o * (P_DERIV_POROSITY * INVERS_FVF_O + P_DERIV_INVERS_FVF_O * POROSITY);
00558
00559 if (FI_CHK_SG (main_vars_, i))
00560 jac_block[p2og_oil_sg] = 0.0;
00561 else if (FI_CHK_RO (main_vars_, i))
00562 {
00563 jac_block[p2og_oil_sg] = VOLUME * POROSITY * local_data_.sat_o * GOR_DERIV_INVERS_FVF;
00564 }
00565 local_data_.ps_block[p2og_oil] = VOLUME * POROSITY * INVERS_FVF_O;
00566 if (FI_CHK_MOMG (main_vars_, i))
00567 {
00568 jac_block[p2og_oil_po] = 0;
00569 local_data_.ps_block[p2og_oil] = VOLUME;
00570 jac_block[p2og_oil_sg] = 0.0;
00571 }
00572 }
00573 else if (n_phases == 2 && is_w)
00574 {
00575 jac_block[p2ow_oil_po] = VOLUME * local_data_.sat_o * (P_DERIV_POROSITY * INVERS_FVF_O + P_DERIV_INVERS_FVF_O * POROSITY);
00576 jac_block[p2ow_oil_so] = VOLUME * POROSITY * INVERS_FVF_O;
00577 local_data_.ps_block[p2ow_oil] = 0.0;
00578 }
00579 else
00580 {
00581 jac_block[0] = VOLUME * (P_DERIV_POROSITY * INVERS_FVF_O + P_DERIV_INVERS_FVF_O * POROSITY);
00582 }
00583 }
00584 }
00585
00590 void
00591 fi_operator_switch_main_vars (double dt)
00592 {
00593 typedef typename calc_model_t::sp_pvt_oil sp_pvt_oil_t;
00594
00595 sp_pvt_oil_t pvt_oil;
00596 index_t prev_pvt_reg = -1;
00597
00598 index_t switch_to_sg_count = 0;
00599 index_t switch_to_ro_count = 0;
00600 index_t switch_to_momg_count = 0;
00601
00602 for (index_t i = 0, cnt = n_cells_; i < cnt; ++i)
00603 {
00604 index_t pvt_reg = pvt_regions_[i];
00605 if (pvt_reg != prev_pvt_reg)
00606 {
00607 prev_pvt_reg = pvt_reg;
00608 pvt_oil = pvt_oil_array [pvt_reg];
00609 }
00610
00611 switch_main_vars <strategy_t>::do_switch (
00612 is_w, is_g, is_o,
00613 d_o, d_g, d_w,
00614 pvt_oil_array [pvt_reg],
00615 pressure_[i],
00616 main_vars_[i],
00617 &saturation_3p_[i * n_phases],
00618 gas_oil_ratio_[i],
00619 calc_model_->lsearch_force_newton_step,
00620 drsdt_,
00621 dt,
00622 calc_model_->old_data_.gas_oil_ratio[i],
00623 switch_to_sg_count,
00624 switch_to_ro_count,
00625 switch_to_momg_count,
00626 i);
00627 }
00628
00629 BOSOUT (section::iters, level::medium)
00630 << "Switch to Sg: " << switch_to_sg_count
00631 << "\tSwitch to Ro: " << switch_to_ro_count
00632 << "\tSwitch to MoMg: " << switch_to_momg_count
00633 << bs_end;
00634 }
00635
00643 BS_FORCE_INLINE void
00644 fill_acc_rhs (index_t i, const data_t &data_i, rhs_item_t *rhs_block, local_data &local_data_)
00645 {
00646
00647 if (is_w)
00648 {
00649 item_t r_wat = 0;
00650 item_t sat_w__ = n_phases > 1 ? local_data_.sat_w : 1.0;
00651
00652 r_wat = -VOLUME * (POROSITY * INVERS_FVF_W * sat_w__ - PREV_FLUID_VOLUME_W);
00653
00654
00655
00656 if (n_phases == 3)
00657 rhs_block[p3_wat] = r_wat;
00658 else if (n_phases == 2)
00659 rhs_block[p2ow_wat] = r_wat;
00660 else
00661 rhs_block[0] = r_wat;
00662 }
00663
00664 if (is_g)
00665 {
00666 item_t r_gas = 0;
00667
00668 if (n_phases > 1)
00669 {
00670 if (FI_CHK_SG (main_vars_, i))
00671 {
00672 r_gas = -VOLUME * (POROSITY * INVERS_FVF_G * local_data_.sat_g + GAS_OIL_RATIO * POROSITY * local_data_.sat_o * INVERS_FVF_O - PREV_FLUID_VOLUME_G);
00673 }
00674 else if (FI_CHK_RO (main_vars_, i))
00675 {
00676 r_gas = -VOLUME * (GAS_OIL_RATIO * POROSITY * local_data_.sat_o * INVERS_FVF_O - PREV_FLUID_VOLUME_G);
00677 }
00678 else if (FI_CHK_MOMG (main_vars_, i))
00679 {
00680 r_gas = VOLUME * PREV_FLUID_VOLUME_G;
00681 }
00682 }
00683 else
00684 {
00685 r_gas = -VOLUME * (POROSITY * INVERS_FVF_G - PREV_FLUID_VOLUME_G);
00686 }
00687
00688
00689
00690 if (n_phases == 3)
00691 rhs_block[p3_gas] = r_gas;
00692 else if (n_phases == 2)
00693 rhs_block[p2og_gas] = r_gas;
00694 else
00695 rhs_block[0] = r_gas;
00696 }
00697
00698 if (is_o)
00699 {
00700 item_t r_oil = 0;
00701 item_t sat_o__ = n_phases > 1 ? local_data_.sat_o : 1.0;
00702
00703 r_oil = -VOLUME * (POROSITY * INVERS_FVF_O * sat_o__ - PREV_FLUID_VOLUME_O);
00704
00705
00706
00707 if (n_phases == 3)
00708 rhs_block[p3_oil] = r_oil;
00709 else if (n_phases == 2 && is_w)
00710 rhs_block[p2ow_oil] = r_oil;
00711 else if (n_phases == 2 && is_g)
00712 rhs_block[p2og_oil] = r_oil;
00713 else
00714 rhs_block[0] = r_oil;
00715 }
00716 }
00717
00726 BS_FORCE_INLINE void
00727 eliminate_cell (const local_data &local_data_, rhs_item_t *jac_block, rhs_item_t *rhs_block, rhs_item_t *sp_block, rhs_item_t *s_rhs_block)
00728 {
00729
00730
00731 M_MINUS_VV_PROD (n_phases, local_data_.ps_block, sp_block, jac_block);
00732
00733 V_MINUS_VS_PROD (n_phases, local_data_.ps_block, s_rhs_block, rhs_block);
00734 }
00735
00739 void
00740 fi_operator_fill ()
00741 {
00742 index_t b_sqr = n_phases * n_phases;
00743 rhs_item_array_t &main_diag_acc = jmatrix_->get_regular_acc_diag();
00744 rhs_item_array_t &ss_diag = jmatrix_->get_ss_diagonal ();
00745 rhs_item_array_t &sp_diag = jmatrix_->get_sp_diagonal ();
00746 rhs_item_array_t &s_rhs = jmatrix_->get_sec_rhs ();
00747 rhs_item_array_t &rhs = jmatrix_->get_rhs ();
00748
00749 BS_ASSERT (!main_diag_acc.empty ());
00750 BS_ASSERT (!ss_diag.empty ());
00751 BS_ASSERT (!sp_diag.empty ());
00752 BS_ASSERT (!s_rhs.empty ());
00753 BS_ASSERT (!rhs.empty ());
00754
00755
00756 for (index_t i = 0; i < n_cells_; ++i)
00757 {
00758 rhs_item_t *jac_block = &main_diag_acc[b_sqr * i];
00759 rhs_item_t *ss_block = &ss_diag[n_sec_vars * n_sec_vars * i];
00760 rhs_item_t *sp_block = &sp_diag[n_sec_vars * n_phases * i];
00761 rhs_item_t *s_rhs_block = &s_rhs[n_sec_vars * i];
00762 rhs_item_t *rhs_block = &rhs[i * n_phases];
00763
00764 index_t i_w = FI_PH_IND (i, d_w, n_phases);
00765 index_t i_g = FI_PH_IND (i, d_g, n_phases);
00766 index_t i_o = FI_PH_IND (i, d_o, n_phases);
00767
00768 local_data local_data_;
00769 const data_t &data_i = data_[i];
00770
00771
00772 if (n_phases > 1)
00773 ss_block[0] = 1;
00774 if (n_phases == 3)
00775 {
00776 s_rhs_block[0] = 1.0 - saturation_3p_[i_w] - saturation_3p_[i_g] - saturation_3p_[i_o];
00777 if (FI_CHK_SG (main_vars_, i))
00778 {
00779 sp_block[p3_sg] = 1.0;
00780 sp_block[p3_so] = 1.0;
00781 sp_block[p3_po] = 0;
00782 }
00783 else if (FI_CHK_RO (main_vars_, i))
00784 {
00785 sp_block[p3_sg] = 0;
00786 sp_block[p3_so] = 1.0;
00787 sp_block[p3_po] = 0;
00788 }
00789 else if (FI_CHK_MOMG (main_vars_, i))
00790 {
00791 sp_block[p3_sg] = 0;
00792 sp_block[p3_so] = 0;
00793 sp_block[p3_po] = 0;
00794 }
00795 }
00796 else if (n_phases == 2 && is_w && is_o)
00797 {
00798 s_rhs_block[0] = 1.0 - saturation_3p_[i_w] - saturation_3p_[i_o];
00799 sp_block[p2ow_so] = 1;
00800 sp_block[p2ow_po] = 0;
00801 }
00802 else if (n_phases == 2 && is_g && is_o)
00803 {
00804 s_rhs_block[0] = 1.0 - saturation_3p_[i_g] - saturation_3p_[i_o];
00805 if (FI_CHK_SG (main_vars_, i))
00806 {
00807 sp_block[p2og_sg] = 1;
00808 sp_block[p2og_po] = 0;
00809 }
00810 else if (FI_CHK_RO (main_vars_, i))
00811 {
00812 sp_block[p2og_sg] = 0;
00813 sp_block[p2og_po] = 0;
00814 }
00815 }
00816
00817 assign (local_data_.ps_block, 0);
00818 local_data_.sat_w = local_data_.sat_g = local_data_.sat_o = 0.0;
00819
00820 if (is_w && n_phases > 1)
00821 local_data_.sat_w = saturation_3p_[i_w];
00822 if (is_g && n_phases > 1)
00823 local_data_.sat_g = saturation_3p_[i_g];
00824 if (is_o && n_phases > 1)
00825 local_data_.sat_o = saturation_3p_[i_o];
00826
00827 fill_jacobian (i, data_i, jac_block, local_data_);
00828 fill_acc_rhs (i, data_i, rhs_block, local_data_);
00829 eliminate_cell (local_data_, jac_block, rhs_block, sp_block, s_rhs_block);
00830
00831
00832
00833
00834
00835
00836
00837
00838
00839
00840
00841
00842
00843
00844 }
00845 }
00846
00862 #ifdef FI_OPERATOR_CELLS_PARALLEL
00863 inline void
00864 fi_operator_cell (index_t istart, const item_t dt,
00865 index_t cell_ind, item_t *jac_block,
00866 item_t *rhs_block,
00867 item_t *ss_block,
00868 item_t *s_rhs_block,
00869 int &switch_to_sg_count,
00870 int &switch_to_ro_count,
00871 int &switch_to_momg_count,
00872 item_t &invers_fvf_average_w,
00873 item_t &invers_fvf_average_g,
00874 item_t &invers_fvf_average_o,
00875 item_t &total_volume,
00876 const sp_mesh_iface_t &mesh)
00877 #else //FI_OPERATOR_CELLS_PARALLEL
00878 inline void
00879 fi_operator_cell (index_t , const item_t dt,
00880 index_t cell_ind, rhs_item_t * ,
00881 rhs_item_t * ,
00882 rhs_item_t * ,
00883 rhs_item_t * ,
00884 int &,
00885 int &,
00886 int &,
00887 item_t &total_volume)
00888 #endif //FI_OPERATOR_CELLS_PARALLEL
00889
00890 {
00891 local_data local_data_;
00892 item_t p_water;
00893 item_t p_gas;
00894 item_t p_oil;
00895 index_t i_temp;
00896
00897 #ifdef _MPI
00898 int n_left = 0;
00899 int n_own = 0;
00900 #endif //_MPI
00901
00902
00903 index_t i = cell_ind;
00904
00905 item_t gor = 0.0;
00906 item_t d_gor = 0.0;
00907
00908 index_t pvt_reg = pvt_regions_[i];
00909
00910 assign (local_data_.ps_block, 0);
00911
00912 if (pressure_[i] < min_p_ || pressure_[i] > max_p_)
00913 {
00914 BOSERR (section::iters, level::error) << boost::format ("Pressure %f in cell [%d] is out of range") %
00915 pressure_[i] % i << bs_end;
00916 }
00917
00918 data_t &data_i = data_[i];
00919
00920
00921 p_oil = p_water = p_gas = pressure_[i];
00922 if (n_phases > 1)
00923 {
00924 if (is_w)
00925 p_water += data_i.cap_pressure[ds_w];
00926 if (is_g)
00927 p_gas += data_i.cap_pressure[ds_g];
00928 }
00929
00930
00931 i_temp = n_phases * i;
00932 if (is_w)
00933 {
00934 pvt_water_array[pvt_reg]->calc(
00935 p_water,
00936 &data_i.invers_fvf[d_w],
00937 &data_i.p_deriv_invers_fvf[d_w],
00938 &data_i.invers_viscosity[d_w],
00939 &data_i.p_deriv_invers_viscosity[d_w],
00940 &data_i.invers_visc_fvf[d_w],
00941 &data_i.p_deriv_invers_visc_fvf[d_w]);
00942 }
00943
00944 if (is_g)
00945 {
00946 pvt_gas_array[pvt_reg]->calc(
00947 p_gas,
00948 &data_i.invers_fvf[d_g],
00949 &data_i.p_deriv_invers_fvf[d_g],
00950 &data_i.invers_viscosity[d_g],
00951 &data_i.p_deriv_invers_viscosity[d_g],
00952 &data_i.invers_visc_fvf[d_g],
00953 &data_i.p_deriv_invers_visc_fvf[d_g]);
00954 }
00955
00956 if (is_o && is_g)
00957 {
00958 pvt_oil_array[pvt_reg]->calc(
00959 is_g, main_vars_[i],p_oil, gas_oil_ratio_[i],
00960 &data_i.invers_fvf[d_o],
00961 &data_i.p_deriv_invers_fvf[d_o],
00962 &data_i.gor_deriv_invers_fvf,
00963 &data_i.invers_viscosity[d_o],
00964 &data_i.p_deriv_invers_viscosity[d_o],
00965 &data_i.gor_deriv_invers_viscosity,
00966 &data_i.invers_visc_fvf[d_o],
00967 &data_i.p_deriv_invers_visc_fvf[d_o],
00968 &data_i.gor_deriv_invers_visc_fvf,
00969 &gor, &d_gor, drsdt_, dt, calc_model_->old_data_.gas_oil_ratio[i]
00970 );
00971 }
00972 else if (is_o)
00973 {
00974 pvt_oil_array[pvt_reg]->calc(
00975 is_g, -1, p_oil, 0,
00976 &data_i.invers_fvf[d_o],
00977 &data_i.p_deriv_invers_fvf[d_o],
00978 0,
00979 &data_i.invers_viscosity[d_o],
00980 &data_i.p_deriv_invers_viscosity[d_o],
00981 0,
00982 &data_i.invers_visc_fvf[d_o],
00983 &data_i.p_deriv_invers_visc_fvf[d_o],
00984 0,
00985 0, 0
00986 );
00987 }
00988
00989
00990 calc_porosity_and_deriv (i,
00991 pvt_reg,
00992 &data_i.porosity,
00993 &data_i.p_deriv_porosity,
00994 &data_i.truns_mult,
00995 &data_i.p_deriv_truns_mult);
00996
00997 total_volume += data_i.porosity * volume_[i];
00998
00999 item_t oil_surface_density = pvt_oil_array[pvt_reg]->get_surface_density ();
01000 item_t gas_surface_density = pvt_gas_array[pvt_reg]->get_surface_density ();
01001 item_t water_surface_density = pvt_water_array[pvt_reg]->get_surface_density ();
01002
01003
01004
01005 if (is_o)
01006 {
01007 if (is_g)
01008 {
01009
01010 if (FI_CHK_SG (main_vars_, i))
01011 {
01012 gas_oil_ratio_[i] = gor;
01013 data_i.p_deriv_gas_oil_ratio = d_gor;
01014 data_i.density[d_o] = data_i.invers_fvf[d_o] * (oil_surface_density + gor * gas_surface_density);
01015
01016 data_i.p_deriv_density[d_o] = data_i.p_deriv_invers_fvf[d_o] * oil_surface_density + gas_surface_density * (gas_oil_ratio_[i] * data_i.p_deriv_invers_fvf[d_o] + data_i.invers_fvf[d_o] * data_i.p_deriv_gas_oil_ratio);
01017 data_i.gor_deriv_density = 0.0;
01018 }
01019 else if (FI_CHK_RO (main_vars_, i))
01020 {
01021 data_i.density[d_o] = data_i.invers_fvf[d_o] * (oil_surface_density + gas_oil_ratio_[i] * gas_surface_density);
01022 data_i.p_deriv_gas_oil_ratio = 0.0;
01023 data_i.p_deriv_density[d_o] = data_i.p_deriv_invers_fvf[d_o] * (oil_surface_density + gas_oil_ratio_[i] * gas_surface_density);
01024 data_i.gor_deriv_density = data_i.gor_deriv_invers_fvf * (oil_surface_density + gas_oil_ratio_[i] * gas_surface_density) + data_i.invers_fvf[d_o] * gas_surface_density;
01025 }
01026 else if (FI_CHK_MOMG (main_vars_, i))
01027 {
01028 gas_oil_ratio_[i] = gor;
01029 data_i.p_deriv_gas_oil_ratio = d_gor;
01030 data_i.density[d_o] = data_i.invers_fvf[d_o] * oil_surface_density;
01031 }
01032 }
01033 else
01034 {
01035 data_i.density[d_o] = data_i.invers_fvf[d_o] * oil_surface_density;
01036 data_i.p_deriv_density[d_o] = data_i.p_deriv_invers_fvf[d_o] * oil_surface_density;
01037 }
01038 }
01039
01040
01041 if (is_w)
01042 {
01043 data_i.density[d_w] = data_i.invers_fvf[d_w] * water_surface_density;
01044 data_i.p_deriv_density[d_w] = data_i.p_deriv_invers_fvf[d_w] * water_surface_density;
01045 }
01046
01047 if (is_g)
01048 {
01049 data_i.density[d_g] = data_i.invers_fvf[d_g] * gas_surface_density;
01050 data_i.p_deriv_density[d_g] = data_i.p_deriv_invers_fvf[d_g] * gas_surface_density;
01051 }
01052
01053
01054
01055
01056 if (is_w)
01057 {
01058
01059 data_i.mobility[d_w] = data_i.invers_visc_fvf[d_w];
01060 data_i.p_deriv_mobility[d_w] = data_i.p_deriv_invers_visc_fvf[d_w];
01061 if (n_phases > 1)
01062 {
01063
01064 data_i.mobility[d_w] *= data_i.relative_perm[d_w];
01065 data_i.p_deriv_mobility[d_w] *= data_i.relative_perm[d_w];
01066 data_i.s_deriv_mobility[d_ww] = data_i.s_deriv_relative_perm[d_ww] * data_i.invers_visc_fvf[d_w]
01067 + data_i.relative_perm[d_w] * data_i.p_deriv_invers_visc_fvf[d_w] * data_i.s_deriv_cap_pressure[ds_w];
01068 if (is_o)
01069 data_i.s_deriv_mobility[d_wo] = 0.0;
01070 if (is_g)
01071 data_i.s_deriv_mobility[d_wg] = 0.0;
01072 }
01073 }
01074
01075 if (is_g)
01076 {
01077 data_i.mobility[d_g] = data_i.invers_visc_fvf[d_g];
01078 data_i.p_deriv_mobility[d_g] = data_i.p_deriv_invers_visc_fvf[d_g];
01079 if (n_phases > 1)
01080 {
01081 data_i.mobility[d_g] *= data_i.relative_perm[d_g];
01082 data_i.p_deriv_mobility[d_g] *= data_i.relative_perm[d_g];
01083
01084 if (FI_CHK_SG (main_vars_, i))
01085 {
01086 data_i.s_deriv_mobility[d_gg] = data_i.s_deriv_relative_perm[d_gg] * data_i.invers_visc_fvf[d_g]
01087 + data_i.relative_perm[d_g] * data_i.p_deriv_invers_visc_fvf[d_g] * data_i.s_deriv_cap_pressure[ds_g];
01088 }
01089 else
01090 {
01091 data_i.s_deriv_mobility[d_gg] = 0;
01092 data_i.mobility[d_g] = 0;
01093 data_i.p_deriv_mobility[d_g] = 0;
01094 }
01095
01096 data_i.s_deriv_mobility[d_go] = 0;
01097
01098 if (is_w)
01099 data_i.s_deriv_mobility[d_gw] = 0;
01100
01101 if (FI_CHK_MOMG (main_vars_, i))
01102 {
01103 data_i.mobility[d_g] = 0;
01104 data_i.p_deriv_mobility[d_g] = 0;
01105 }
01106 }
01107 }
01108
01109 if (is_o)
01110 {
01111 data_i.mobility[d_o] = data_i.invers_visc_fvf[d_o];
01112 data_i.p_deriv_mobility[d_o] = data_i.p_deriv_invers_visc_fvf[d_o];
01113 if (n_phases > 1)
01114 {
01115 data_i.mobility[d_o] *= data_i.relative_perm[d_o];
01116 data_i.p_deriv_mobility[d_o] *= data_i.relative_perm[d_o];
01117
01118 data_i.s_deriv_mobility[d_oo] = data_i.s_deriv_relative_perm[d_oo] * data_i.invers_visc_fvf[d_o];
01119
01120 if (is_w)
01121 data_i.s_deriv_mobility[d_ow] = data_i.s_deriv_relative_perm[d_ow] * data_i.invers_visc_fvf[d_o];
01122
01123 if (is_g)
01124 {
01125 if (FI_CHK_SG (main_vars_, i))
01126 data_i.s_deriv_mobility[d_og] = data_i.s_deriv_relative_perm[d_og] * data_i.invers_visc_fvf[d_o];
01127 else if (FI_CHK_RO (main_vars_, i))
01128 data_i.s_deriv_mobility[d_og] = data_i.relative_perm[d_o] * data_i.gor_deriv_invers_visc_fvf;
01129 else if (FI_CHK_MOMG (main_vars_, i))
01130 {
01131 data_i.mobility[d_o] = 0;
01132 data_i.p_deriv_mobility[d_o] = 0;
01133 data_i.s_deriv_mobility[d_oo] = 0;
01134 data_i.s_deriv_mobility[d_og] = 0;
01135 if (is_w)
01136 data_i.s_deriv_mobility[d_ow] = 0;
01137 }
01138 }
01139 }
01140 }
01141
01142 #ifdef FI_OPERATOR_CELLS_PARALLEL
01143 if (is_w)
01144 invers_fvf_average_w += invers_fvf[i_w];
01145 if (is_g)
01146 invers_fvf_average_g += invers_fvf[i_g];
01147 if (is_o)
01148 invers_fvf_average_o += invers_fvf[i_o];
01149 #else //FI_OPERATOR_CELLS_PARALLEL
01150
01151 #ifdef _MPI
01152 if (i >= n_left && i < n_own)
01153 {
01154 #endif //_MPI
01155 if (is_w)
01156 calc_model_->invers_fvf_average[d_w] += data_i.invers_fvf[d_w];
01157 if (is_g)
01158 calc_model_->invers_fvf_average[d_g] += data_i.invers_fvf[d_g];
01159 if (is_o)
01160 calc_model_->invers_fvf_average[d_o] += data_i.invers_fvf[d_o];
01161 #ifdef _MPI
01162 }
01163 #endif //_MPI
01164
01165 #endif //FI_OPERATOR_CELLS_PARALLEL
01166
01167 }
01168
01174 bool
01175 check_norm (index_t istart)
01176 {
01177 if (!istart)
01178 {
01179 #ifdef _DEBUG
01180 BOSOUT (section::iters, level::debug) << "rhs_residual: " << fabs (norm_.val[norms::C_CPV]) << " / " << rhs_residual_ << bs_end;
01181 BOSOUT (section::iters, level::debug) << "mb_error: " << fabs (norm_.val[norms::MB_ERR]) << " / " << mb_error_ << bs_end;
01182 BOSOUT (section::iters, level::debug) << "s_rhs: " << fabs (norm_.val[norms::S_RHS]) << " / " << s_rhs_norm_ << bs_end;
01183 #endif
01184
01185 return fabs (norm_.val[norms::C_CPV]) < rhs_residual_
01186 && fabs (norm_.val[norms::MB_ERR]) < mb_error_
01187
01188 ;
01189 }
01190
01191 return false;
01192 }
01197 bool
01198 check_solution_mult_cell (index_t istart, index_t iteration, index_t n_approx, const norms_storage_t &base_norm)
01199 {
01200
01201
01202 return false;
01203
01204 if (!istart && !iteration)
01205 {
01206 BOSOUT (section::iters, level::debug) << "check_solution_mult_cell" << bs_end;
01207
01208 item_t mult = calc_solution_mult_cell (base_norm);
01209 if (mult < 0.9 && iteration < n_approx - 1)
01210 {
01211 if (mult < item_t (0.1))
01212 mult = item_t (0.1);
01213
01214 BOSOUT (section::iters, level::debug) << "check_solution_mult_cell 2, mult = " << mult << bs_end;
01215
01216 restore_prev_niter_vars ();
01217 if (calc_model_->apply_newton_correction (mult, 2, mesh_, jmatrix_))
01218 bs_throw_exception ("apply_newton_correction failed");
01219
01220 return true;
01221 }
01222 #ifdef _DEBUG
01223 else
01224 {
01225 BOSOUT (section::iters, level::debug) << "check_solution_mult_cell 3, mult = " << mult << bs_end;
01226 }
01227 #endif
01228 }
01229
01230 return false;
01231 }
01232
01238 void
01239 do_dt_reduce (double &dt, index_t istart)
01240 {
01241 if (istart)
01242 {
01243 item_t dt_mult = calc_step_dt_mult (0, (item_t) calc_model_->ts_params->get_float (fi_params::MAX_NORM_ON_FIRST_N));
01244 if (dt_mult < 0.95)
01245 {
01246 jmatrix_->mult_flux_part (dt_mult);
01247 dt *= dt_mult;
01248
01249 BOSOUT (section::iters, level::medium) << "dT reduced by factor " << dt_mult << ": " << dt << bs_end;
01250 }
01251 dt_mult = 0;
01252 }
01253 }
01259 void
01260 do_dt_tuning (const double &dt, bool tuning)
01261 {
01262 if (tuning)
01263 {
01264 item_t dt_mult = calc_step_dt_mult (0, (item_t) calc_model_->ts_params->get_float (fi_params::MAX_NORM_ON_TS));
01265 if (dt_mult < 0.95)
01266 {
01267 jmatrix_->mult_flux_part (dt_mult);
01268
01269 BOSOUT (section::iters, level::medium) << "Use dT " << dt_mult * (dt) << " instead of original dT " << dt << bs_end;
01270 }
01271 else
01272 dt_mult = 1.0;
01273 }
01274 }
01275
01279 void
01280 save_prev_niter_vars ()
01281 {
01282 int is_line_search = calc_model_->ts_params->get_int (fi_params::SELECT_SOL_STEPS);
01283 if (is_line_search != 0)
01284 {
01285 calc_model_->prev_niter_data_.save (calc_model_);
01286 reservoir_->pre_newton_step ();
01287 }
01288 }
01289
01293 void
01294 restore_prev_niter_vars ()
01295 {
01296 reservoir_->restart_newton_step ();
01297 }
01298
01302 void
01303 norm_calc ();
01304
01310 void
01311 update_norm_by_cell (index_t i, norms_storage_t &ns);
01312
01319 item_t
01320 calc_step_dt_mult (item_t prev_mult, item_t max_res);
01321
01333 void
01334 calc_porosity_and_deriv (index_t i,
01335 index_t pvt_reg,
01336 item_t *poro,
01337 item_t *dp_poro,
01338 item_t *t_mult,
01339 item_t *dp_t_mult);
01340
01344 void
01345 calc_prev_fluid_volume ();
01346
01351 void
01352 debug_save_data (item_t dt);
01353
01359 item_t
01360 calc_solution_mult_cell (const norms_storage_t &base_norm);
01361
01367 bool
01368 block_connections_mpfa (const item_t &dt);
01369
01370 public:
01371 sp_calc_model_t &calc_model_;
01372 const sp_rock_grid_prop_t &rock_grid_prop_;
01373 sp_reservoir_t &reservoir_;
01374 sp_jacobian_t &jacobian_;
01375 sp_jmatrix_t &jmatrix_;
01376 const sp_mesh_iface_t &mesh_;
01377
01378 sp_bcsr_matrix_t trns_matrix_;
01379 const rhs_item_array_t &trns_values_;
01380 const index_array_t &trns_rows_ptr_;
01381 const index_array_t &trns_cols_ptr_;
01382
01383 sp_bcsr_matrix_t reg_matrix_;
01384 rhs_item_array_t ®_values_;
01385 const index_array_t ®_rows_ptr_;
01386 const index_array_t ®_cols_ptr_;
01387
01388 const index_array_t &m_array_;
01389 const index_array_t &p_array_;
01390
01391 rhs_item_array_t &rhs_;
01392 item_array_t &sol_;
01393 rhs_item_array_t &flux_rhs_;
01394 rhs_item_array_t &sp_diag_;
01395 rhs_item_array_t &s_rhs_;
01396
01397 const item_array_t &depths_;
01398
01399 index_t n_cells_;
01400 index_t n_connections_;
01401
01402 index_t n_sec_vars;
01403
01404 index_t d_w;
01405 index_t d_g;
01406 index_t d_o;
01407 index_t ds_w;
01408 index_t ds_g;
01409
01410 index_t d_gg, d_gw, d_go;
01411 index_t d_wg, d_ww, d_wo;
01412 index_t d_og, d_ow, d_oo;
01413
01414 data_array_t &data_;
01415 item_array_t &saturation_3p_;
01416 const item_array_t &pressure_;
01417 item_array_t &gas_oil_ratio_;
01418 main_var_array_t &main_vars_;
01419 const item_array_t &volume_;
01420 const item_array_t &poro_array_;
01421 const item_array_t &rock_grid_comp_const_;
01422 const item_array_t &rock_grid_comp_ref_pressure_;
01423 const index_array_t &sat_regions_;
01424 const index_array_t &pvt_regions_;
01425
01426 const sp_pvt_dead_oil_array_t &pvt_oil_array;
01427 const sp_pvt_water_array_t &pvt_water_array;
01428 const sp_pvt_gas_array_t &pvt_gas_array;
01429
01430 item_t min_p_;
01431 item_t max_p_;
01432 item_t drsdt_;
01433 item_t rhs_residual_;
01434 item_t mb_error_;
01435 item_t s_rhs_norm_;
01436
01437 norms_storage_t &norm_;
01438 rhs_item_array_t &cfl_;
01439 };
01440
01441 }
01442
01443 #include "fi_operator_norm_calc.h"
01444 #include "fi_operator_calc_step_dt_mult.h"
01445 #include "fi_operator_calc_porosity_and_deriv.h"
01446 #include "fi_operator_calc_prev_fluid_volume.h"
01447 #include "fi_operator_save_debug_data.h"
01448 #include "fi_operator_calc_solution_mult_cell.h"
01449
01450 #ifdef BS_USE_TPFA_
01451 #include "fi_operator_block_connections_mpfa.h"
01452 #else
01453 #include "fi_operator_block_connections_mpfa_2.h"
01454 #endif
01455
01456 #endif // #ifndef BS_FI_OPERATOR_H_
01457
