00001
00009 #ifndef BS_FI_OPERATOR_BLOCK_CONNECTION_MPFA_2_H_
00010 #define BS_FI_OPERATOR_BLOCK_CONNECTION_MPFA_2_H_
00011
00012 #include "calc_model.h"
00013 #include "matrix_vector_op.h"
00014 #include "pp_index.h"
00015
00016 #define IS_STABLE 0
00017 #define UNKNOWN -1
00018
00019 namespace blue_sky {
00020
00021 namespace mpfa
00022 {
00023 #define PSI_W psi_[d_w]
00024 #define PSI_G psi_[d_g]
00025 #define PSI_O psi_[d_o]
00026 #define RHO_W rho_[d_w]
00027 #define RHO_G rho_[d_g]
00028 #define RHO_O rho_[d_o]
00029 #define UP_W up_[d_w]
00030 #define UP_G up_[d_g]
00031 #define UP_O up_[d_o]
00032 #define UP_CELL_W up_cell_[d_w]
00033 #define UP_CELL_G up_cell_[d_g]
00034 #define UP_CELL_O up_cell_[d_o]
00035 #define FLOW_W flow_[d_w]
00036 #define FLOW_G flow_[d_g]
00037 #define FLOW_O flow_[d_o]
00038
00039 #define CAP_PRESSURE_K_W data_[k_cell].cap_pressure [d_w]
00040 #define CAP_PRESSURE_K_G data_[k_cell].cap_pressure [d_g]
00041 #define PRESSURE_K pressure_[k_cell]
00042 #define S_DERIV_CAP_PRESSURE_K_W data_[k_cell].s_deriv_cap_pressure [ds_w]
00043 #define S_DERIV_CAP_PRESSURE_K_G data_[k_cell].s_deriv_cap_pressure [ds_g]
00044
00045 #define MOBILITY_UP_W mobility_up_w
00046 #define MOBILITY_UP_G mobility_up_g
00047 #define MOBILITY_UP_O mobility_up_o
00048 #define GOR_UP_O gor_up_o
00049 #define SW_DERIV_MOBILITY_UP_W sw_deriv_mobility_up_w
00050 #define SW_DERIV_MOBILITY_UP_G sw_deriv_mobility_up_g
00051 #define SW_DERIV_MOBILITY_UP_O sw_deriv_mobility_up_o
00052 #define SG_DERIV_MOBILITY_UP_W sg_deriv_mobility_up_w
00053 #define SG_DERIV_MOBILITY_UP_G sg_deriv_mobility_up_g
00054 #define SG_DERIV_MOBILITY_UP_O sg_deriv_mobility_up_o
00055 #define SO_DERIV_MOBILITY_UP_W so_deriv_mobility_up_w
00056 #define SO_DERIV_MOBILITY_UP_G so_deriv_mobility_up_g
00057 #define SO_DERIV_MOBILITY_UP_O so_deriv_mobility_up_o
00058 #define P_DERIV_MOBILITY_UP_W p_deriv_mobility_up_w
00059 #define P_DERIV_MOBILITY_UP_G p_deriv_mobility_up_g
00060 #define P_DERIV_MOBILITY_UP_O p_deriv_mobility_up_o
00061 #define MAIN_VAR_UP_G main_var_up_g
00062 #define MAIN_VAR_UP_O main_var_up_o
00063 #define P_DERIV_GOR_UP_O p_deriv_gor_up_o
00064
00073 template <typename strategy_t, bool is_w, bool is_g, bool is_o>
00074 struct mpfa_base_impl
00075 {
00076 public:
00077 typedef typename strategy_t::item_t item_t;
00078 typedef typename strategy_t::rhs_item_t rhs_item_t;
00079 typedef typename strategy_t::index_t index_t;
00080 typedef typename strategy_t::item_array_t item_array_t;
00081 typedef typename strategy_t::rhs_item_array_t rhs_item_array_t;
00082 typedef typename strategy_t::index_array_t index_array_t;
00083 typedef boost::array <index_t, FI_PHASE_TOT> up_cell_array_t;
00084 typedef calc_model <strategy_t> calc_model_t;
00085 typedef typename calc_model_t::data_t data_t;
00086 typedef typename calc_model_t::data_array_t data_array_t;
00087 typedef typename calc_model_t::main_var_array_t main_var_array_t;
00088
00089 enum
00090 {
00091 n_phases = is_w + is_g + is_o,
00092 b_sqr = n_phases * n_phases,
00093 };
00094
00095 enum
00096 {
00097 is_1p = n_phases == 1,
00098 is_2p = n_phases == 2,
00099 is_3p = n_phases == 3,
00100 };
00101
00102 enum
00103 {
00104 gas_sg = detail::pp_index <n_phases, is_w, is_g>::gas_sg,
00105 gas_so = detail::pp_index <n_phases, is_w, is_g>::gas_so,
00106 gas_po = detail::pp_index <n_phases, is_w, is_g>::gas_po,
00107 oil_sg = detail::pp_index <n_phases, is_w, is_g>::oil_sg,
00108 oil_so = detail::pp_index <n_phases, is_w, is_g>::oil_so,
00109 oil_po = detail::pp_index <n_phases, is_w, is_g>::oil_po,
00110 wat_sg = detail::pp_index <n_phases, is_w, is_g>::wat_sg,
00111 wat_so = detail::pp_index <n_phases, is_w, is_g>::wat_so,
00112 wat_po = detail::pp_index <n_phases, is_w, is_g>::wat_po,
00113 };
00114
00115 enum
00116 {
00117 gas_idx = 0,
00118 oil_idx = is_g,
00119 wat_idx = is_g + is_o,
00120
00121 gas_sw = !is_1p ? b_sqr + gas_idx : -1,
00122 oil_sw = !is_1p ? b_sqr + oil_idx : -1,
00123 wat_sw = !is_1p ? b_sqr + wat_idx : -1,
00124 };
00125
00126 public:
00133 mpfa_base_impl (const fi_operator_impl <strategy_t, is_w, is_g, is_o> &fi_operator_, rhs_item_array_t &rhs_, rhs_item_array_t ®_values_)
00134 : d_w (fi_operator_.d_w),
00135 d_g (fi_operator_.d_g),
00136 d_o (fi_operator_.d_o),
00137 ds_w (fi_operator_.ds_w),
00138 ds_g (fi_operator_.ds_g),
00139 saturation_3p_ (fi_operator_.saturation_3p_),
00140 pressure_ (fi_operator_.pressure_),
00141 gas_oil_ratio_ (fi_operator_.gas_oil_ratio_),
00142 main_vars_ (fi_operator_.main_vars_),
00143 data_ (fi_operator_.data_),
00144 depths_ (fi_operator_.depths_),
00145 gravity_ (fi_operator_.calc_model_->internal_constants.gravity_constant),
00146 n_sec_vars (fi_operator_.n_sec_vars),
00147 sp_diag_ (fi_operator_.sp_diag_),
00148 s_rhs_ (fi_operator_.s_rhs_),
00149 reg_values_ (reg_values_),
00150 rhs_ (rhs_),
00151 m_mem (0),
00152 p_mem (0),
00153 depth_top (0)
00154 {
00155 }
00156
00161 BS_FORCE_INLINE void
00162 fill_jacobian_k_derivs (index_t k) const
00163 {
00164 rhs_item_t *jacobian_ik = ®_values_ [m_mem[k] * b_sqr];
00165 rhs_item_t *jacobian_jk = ®_values_ [p_mem[k] * b_sqr];
00166
00167 for (index_t i = 0; i < b_sqr; ++i)
00168 {
00169 jacobian_ik [i] += pp[i];
00170 jacobian_jk [i] -= pp[i];
00171 }
00172 }
00173
00179 void
00180 compute_rho (const index_t *cells, index_t cells_count)
00181 {
00182 boost::array <item_t, n_phases> counter;
00183 boost::array <index_t, n_phases> flag;
00184
00185 assign (rho_, 0);
00186 assign (counter, 0);
00187 assign (flag, 0);
00188
00189 for (index_t j = 0; j < cells_count; ++j)
00190 {
00191 index_t ic = cells[j];
00192 index_t idx = ic * n_phases;
00193 const data_t &data_ic = data_[ic];
00194
00195 if (n_phases > 1)
00196 {
00197 bool is_g_ = is_g && saturation_3p_[idx + d_g] > EPS_DIV;
00198 bool is_o_ = is_o && saturation_3p_[idx + d_o] > EPS_DIV;
00199 bool is_w_ = is_w && saturation_3p_[idx + d_w] > EPS_DIV;
00200
00201 if (is_g_)
00202 {
00203 rho_[d_g] += data_ic.density[d_g];
00204 ++counter[d_g];
00205 ++flag[d_g];
00206 }
00207 if (is_o_)
00208 {
00209 rho_[d_o] += data_ic.density[d_o];
00210 ++counter[d_o];
00211 ++flag[d_o];
00212 }
00213 if (is_w_)
00214 {
00215 rho_[d_w] += data_ic.density[d_w];
00216 ++counter[d_w];
00217 ++flag[d_w];
00218 }
00219 }
00220 else
00221 {
00222 rho_[0] += data_ic.density[0];
00223 ++counter[0];
00224 ++flag[0];
00225 }
00226 }
00227 for (index_t j = 0; j < n_phases; ++j)
00228 {
00229 if (flag[j])
00230 rho_[j] /= counter[j];
00231 }
00232 }
00233
00240 void
00241 compute_psi (const rhs_item_t *truns, const index_t *cells, index_t cells_count)
00242 {
00243 assign (psi_, 0);
00244
00245 depth_top = depths_ [cells[0]];
00246 for (index_t k = 0; k < cells_count; ++k)
00247 {
00248 index_t k_cell = cells [k];
00249 rhs_item_t truns_k = truns[k];
00250 item_t g_h = gravity_ * (depths_[k_cell] - depth_top);
00251 item_t p = PRESSURE_K;
00252
00253 if (is_g) PSI_G += truns_k * (p + CAP_PRESSURE_K_G - RHO_G * g_h);
00254 if (is_o) PSI_O += truns_k * (p - RHO_O * g_h);
00255 if (is_w) PSI_W += truns_k * (p + CAP_PRESSURE_K_W - RHO_W * g_h);
00256 }
00257 }
00258
00267 void
00268 compute_upstream (index_t i, index_t j, index_t i_cell, index_t j_cell)
00269 {
00270 index_t up_i[] = {j, i};
00271 index_t up_cell[] = {j_cell, i_cell};
00272
00273 if (is_g) UP_G = up_i[PSI_G > 0.0];
00274 if (is_o) UP_O = up_i[PSI_O > 0.0];
00275 if (is_w) UP_W = up_i[PSI_W > 0.0];
00276 if (is_g) UP_CELL_G = up_cell[PSI_G > 0.0];
00277 if (is_o) UP_CELL_O = up_cell[PSI_O > 0.0];
00278 if (is_w) UP_CELL_W = up_cell[PSI_W > 0.0];
00279
00280 if (is_w) sw_deriv_mobility_up_w = data_[UP_CELL_W].s_deriv_mobility[d_w * n_phases + d_w];
00281 if (is_w && is_g) sw_deriv_mobility_up_g = data_[UP_CELL_G].s_deriv_mobility[d_w * n_phases + d_g];
00282 if (is_w && is_o) sw_deriv_mobility_up_o = data_[UP_CELL_O].s_deriv_mobility[d_w * n_phases + d_o];
00283 if (is_g && is_w) sg_deriv_mobility_up_w = data_[UP_CELL_W].s_deriv_mobility[d_g * n_phases + d_w];
00284 if (is_g) sg_deriv_mobility_up_g = data_[UP_CELL_G].s_deriv_mobility[d_g * n_phases + d_g];
00285 if (is_g && is_o) sg_deriv_mobility_up_o = data_[UP_CELL_O].s_deriv_mobility[d_g * n_phases + d_o];
00286 if (is_o && is_w) so_deriv_mobility_up_w = data_[UP_CELL_W].s_deriv_mobility[d_o * n_phases + d_w];
00287 if (is_o && is_g) so_deriv_mobility_up_g = data_[UP_CELL_G].s_deriv_mobility[d_o * n_phases + d_g];
00288 if (is_o) so_deriv_mobility_up_o = data_[UP_CELL_O].s_deriv_mobility[d_o * n_phases + d_o];
00289 if (is_w) p_deriv_mobility_up_w = data_[UP_CELL_W].p_deriv_mobility[d_w];
00290 if (is_g) p_deriv_mobility_up_g = data_[UP_CELL_G].p_deriv_mobility[d_g];
00291 if (is_o) p_deriv_mobility_up_o = data_[UP_CELL_O].p_deriv_mobility[d_o];
00292 if (is_w) mobility_up_w = data_[UP_CELL_W].mobility [d_w];
00293 if (is_g) mobility_up_g = data_[UP_CELL_G].mobility [d_g];
00294 if (is_o) mobility_up_o = data_[UP_CELL_O].mobility [d_o];
00295 if (is_o && is_g) gor_up_o = gas_oil_ratio_[UP_CELL_O];
00296 if (is_o && is_g) p_deriv_gor_up_o = data_[UP_CELL_O].p_deriv_gas_oil_ratio;
00297 if (is_g) main_var_up_g = main_vars_[UP_CELL_G];
00298 if (is_o) main_var_up_o = main_vars_[UP_CELL_O];
00299 }
00300
00305 void
00306 compute_flow (const double &dt)
00307 {
00308 if (is_g) FLOW_G = dt * MOBILITY_UP_G * PSI_G;
00309 if (is_o) FLOW_O = dt * MOBILITY_UP_O * PSI_O;
00310 if (is_w) FLOW_W = dt * MOBILITY_UP_W * PSI_W;
00311
00312 if (is_g && is_o)
00313 FLOW_G += FLOW_O * GOR_UP_O;
00314 }
00315
00321 void
00322 fill_rhs (index_t i_cell, index_t j_cell)
00323 {
00324 if (is_g)
00325 {
00326 rhs_ [i_cell * n_phases + gas_idx] += -FLOW_G;
00327 rhs_ [j_cell * n_phases + gas_idx] -= -FLOW_G;
00328 }
00329 if (is_o)
00330 {
00331 rhs_ [i_cell * n_phases + oil_idx] += -FLOW_O;
00332 rhs_ [j_cell * n_phases + oil_idx] -= -FLOW_O;
00333 }
00334 if (is_w)
00335 {
00336 rhs_ [i_cell * n_phases + wat_idx] += -FLOW_W;
00337 rhs_ [j_cell * n_phases + wat_idx] -= -FLOW_W;
00338 }
00339 }
00340
00348 inline item_t
00349 wat_sw_deriv (index_t k_cell, main_var_type main_var, item_t truns)
00350 {
00351 return MOBILITY_UP_W * truns * S_DERIV_CAP_PRESSURE_K_W;
00352 }
00357 inline item_t
00358 wat_sw_deriv_up ()
00359 {
00360 return SW_DERIV_MOBILITY_UP_W * PSI_W;
00361 }
00369 inline item_t
00370 wat_po_deriv (index_t k_cell, main_var_type main_var, item_t truns)
00371 {
00372 return MOBILITY_UP_W * truns;
00373 }
00378 inline item_t
00379 wat_po_deriv_up ()
00380 {
00381 return P_DERIV_MOBILITY_UP_W * PSI_W;
00382 }
00387 inline item_t
00388 wat_sg_deriv_up ()
00389 {
00390 return SW_DERIV_MOBILITY_UP_G * PSI_W;
00391 }
00396 inline item_t
00397 wat_so_deriv_up ()
00398 {
00399 return SW_DERIV_MOBILITY_UP_O * PSI_W;
00400 }
00401
00406 inline item_t
00407 gas_sw_deriv_up ()
00408 {
00409 return SG_DERIV_MOBILITY_UP_W * PSI_G;
00410 }
00418 inline item_t
00419 gas_po_deriv (index_t k_cell, main_var_type main_var, item_t truns)
00420 {
00421 return MOBILITY_UP_G * truns;
00422 }
00427 inline item_t
00428 gas_po_deriv_up ()
00429 {
00430 return P_DERIV_MOBILITY_UP_G * PSI_G;
00431 }
00439 inline item_t
00440 gas_sg_deriv (index_t k_cell, main_var_type main_var, item_t truns)
00441 {
00442 return main_var == FI_SG_VAR ? MOBILITY_UP_G * truns * S_DERIV_CAP_PRESSURE_K_G : 0;
00443 }
00448 inline item_t
00449 gas_sg_deriv_up ()
00450 {
00451 return MAIN_VAR_UP_G == FI_SG_VAR ? SG_DERIV_MOBILITY_UP_G * PSI_G : 0;
00452 }
00457 inline item_t
00458 gas_so_deriv_up ()
00459 {
00460 return SG_DERIV_MOBILITY_UP_O * PSI_G;
00461 }
00462
00467 inline item_t
00468 oil_sw_deriv_up ()
00469 {
00470 return SO_DERIV_MOBILITY_UP_W * PSI_O;
00471 }
00479 inline item_t
00480 oil_po_deriv (index_t k_cell, main_var_type main_var, item_t truns)
00481 {
00482 return MOBILITY_UP_O * truns;
00483 }
00488 inline item_t
00489 oil_po_deriv_up ()
00490 {
00491 return P_DERIV_MOBILITY_UP_O * PSI_O;
00492 }
00497 inline item_t
00498 oil_sg_deriv_up ()
00499 {
00500 return MAIN_VAR_UP_O == FI_SG_VAR ? (SO_DERIV_MOBILITY_UP_G * PSI_O) : (data_[UP_CELL_O].gor_deriv_invers_visc_fvf * data_[UP_CELL_O].relative_perm[d_o] * PSI_O);
00501 }
00506 inline item_t
00507 oil_so_deriv_up ()
00508 {
00509 return SO_DERIV_MOBILITY_UP_O * PSI_O;
00510 }
00511
00542 void
00543 compute_k_derivs (index_t k, index_t k_cell, main_var_type main_var, item_t truns, index_t i_cell, index_t j_cell)
00544 {
00545 assign (pp, 0);
00546
00547
00548 if (is_g && !is_1p) pp [gas_sg] = gas_sg_deriv (k_cell, main_var, truns);
00549 if (is_g) pp [gas_po] = gas_po_deriv (k_cell, main_var, truns);
00550 if (is_o) pp [oil_po] = oil_po_deriv (k_cell, main_var, truns);
00551 if (is_w) pp [wat_po] = wat_po_deriv (k_cell, main_var, truns);
00552 if (is_w && !is_1p) pp [wat_sw] = wat_sw_deriv (k_cell, main_var, truns);
00553
00554 if (is_g && is_o) pp [gas_po] += GOR_UP_O * pp [oil_po];
00555
00556 if (!is_1p)
00557 {
00558 const rhs_item_t *sp = &sp_diag_ [k_cell * n_phases];
00559 const item_t sr = s_rhs_ [k_cell];
00560 rhs_item_t *ri = &rhs_ [i_cell * n_phases];
00561 rhs_item_t *rj = &rhs_ [j_cell * n_phases];
00562
00563 m_minus_vv_prod <n_phases>::eliminate (&pp[b_sqr], sp, pp);
00564 v_minus_vs_prod <n_phases>::eliminate (&pp[b_sqr], sr, ri);
00565 v_minus_vs_prod <n_phases>::eliminate (&pp[b_sqr], -sr, rj);
00566 }
00567 fill_jacobian_k_derivs (k);
00568 }
00569
00577 void
00578 compute_up_derivs (const double &dt, index_t i_cell, index_t j_cell)
00579 {
00580 if (is_g && !is_1p) pp [gas_sg] = dt * gas_sg_deriv_up ();
00581 if (is_g && !is_1p) pp [gas_so] = dt * gas_so_deriv_up ();
00582 if (is_g) pp [gas_po] = dt * gas_po_deriv_up ();
00583
00584 if (is_o && is_g) pp [oil_sg] = dt * oil_sg_deriv_up ();
00585 if (is_o && !is_1p) pp [oil_so] = dt * oil_so_deriv_up ();
00586 if (is_o) pp [oil_po] = dt * oil_po_deriv_up ();
00587
00588 if (is_w && is_g) pp [wat_sg] = dt * wat_sg_deriv_up ();
00589 if (is_w && is_o) pp [wat_so] = dt * wat_so_deriv_up ();
00590 if (is_w) pp [wat_po] = dt * wat_po_deriv_up ();
00591
00592 if (is_g && is_w) pp [gas_sw] = dt * gas_sw_deriv_up ();
00593 if (is_o && is_w) pp [oil_sw] = dt * oil_sw_deriv_up ();
00594 if (is_w && !is_1p) pp [wat_sw] = dt * wat_sw_deriv_up ();
00595
00596 if (!is_1p)
00597 eliminate (i_cell, j_cell);
00598
00599 fill_jacobian ();
00600
00601 if (is_g && is_o)
00602 {
00603 pp [gas_sg] = GOR_UP_O * dt * oil_sg_deriv_up ();
00604 pp [gas_so] = GOR_UP_O * dt * oil_so_deriv_up ();
00605 pp [gas_po] = GOR_UP_O * dt * oil_po_deriv_up ();
00606 pp [gas_sw] = GOR_UP_O * pp [oil_sw];
00607
00608 if (MAIN_VAR_UP_O == FI_SG_VAR)
00609 {
00610 pp [gas_po] += P_DERIV_GOR_UP_O * FLOW_O;
00611 }
00612 else
00613 {
00614 pp [gas_sg] += FLOW_O;
00615 }
00616
00617 BS_ASSERT (gas_sg == 0) ((int)gas_sg);
00618 BS_ASSERT (p3_gas == 0) ((int)p3_gas);
00619
00620 if (!is_1p)
00621 {
00622 const rhs_item_t *sp_o = &sp_diag_ [UP_CELL_O * n_phases];
00623 pp [gas_idx * n_phases + 0] -= pp [gas_sw] * sp_o [0];
00624 pp [gas_idx * n_phases + 1] -= pp [gas_sw] * sp_o [1];
00625 pp [gas_idx * n_phases + 2] -= pp [gas_sw] * sp_o [2];
00626
00627 rhs_ [i_cell * n_phases + gas_idx] -= pp [gas_sw] * s_rhs_ [UP_CELL_O];
00628 rhs_ [j_cell * n_phases + gas_idx] += pp [gas_sw] * s_rhs_ [UP_CELL_O];
00629 }
00630
00631 for (index_t i = 0; i < n_phases; ++i)
00632 {
00633 reg_values_ [i + m_mem[UP_O] * b_sqr + gas_idx * n_phases] += pp[i + gas_idx * n_phases];
00634 reg_values_ [i + p_mem[UP_O] * b_sqr + gas_idx * n_phases] -= pp[i + gas_idx * n_phases];
00635 }
00636 }
00637 }
00638
00644 void
00645 eliminate (index_t i_cell, index_t j_cell)
00646 {
00647 BS_ASSERT (!is_1p);
00648 if (is_g)
00649 {
00650 const rhs_item_t *sp_g = &sp_diag_ [UP_CELL_G * n_phases];
00651 if (n_phases > 0) pp [gas_idx * n_phases + 0] -= pp [gas_sw] * sp_g [0];
00652 if (n_phases > 1) pp [gas_idx * n_phases + 1] -= pp [gas_sw] * sp_g [1];
00653 if (n_phases > 2) pp [gas_idx * n_phases + 2] -= pp [gas_sw] * sp_g [2];
00654 }
00655
00656 if (is_o)
00657 {
00658 const rhs_item_t *sp_o = &sp_diag_ [UP_CELL_O * n_phases];
00659 if (n_phases > 0) pp [oil_idx * n_phases + 0] -= pp [oil_sw] * sp_o [0];
00660 if (n_phases > 1) pp [oil_idx * n_phases + 1] -= pp [oil_sw] * sp_o [1];
00661 if (n_phases > 2) pp [oil_idx * n_phases + 2] -= pp [oil_sw] * sp_o [2];
00662 }
00663
00664 if (is_w)
00665 {
00666 const rhs_item_t *sp_w = &sp_diag_ [UP_CELL_W * n_phases];
00667 if (n_phases > 0) pp [wat_idx * n_phases + 0] -= pp [wat_sw] * sp_w [0];
00668 if (n_phases > 1) pp [wat_idx * n_phases + 1] -= pp [wat_sw] * sp_w [1];
00669 if (n_phases > 2) pp [wat_idx * n_phases + 2] -= pp [wat_sw] * sp_w [2];
00670 }
00671
00672 if (is_g) rhs_ [i_cell * n_phases + gas_idx] -= pp [gas_sw] * s_rhs_ [UP_CELL_G];
00673 if (is_o) rhs_ [i_cell * n_phases + oil_idx] -= pp [oil_sw] * s_rhs_ [UP_CELL_O];
00674 if (is_w) rhs_ [i_cell * n_phases + wat_idx] -= pp [wat_sw] * s_rhs_ [UP_CELL_W];
00675
00676 if (is_g) rhs_ [j_cell * n_phases + gas_idx] += pp [gas_sw] * s_rhs_ [UP_CELL_G];
00677 if (is_o) rhs_ [j_cell * n_phases + oil_idx] += pp [oil_sw] * s_rhs_ [UP_CELL_O];
00678 if (is_w) rhs_ [j_cell * n_phases + wat_idx] += pp [wat_sw] * s_rhs_ [UP_CELL_W];
00679 }
00680
00684 void
00685 fill_jacobian ()
00686 {
00687 for (index_t i = 0; i < n_phases; ++i)
00688 {
00689 if (is_g) reg_values_ [i + m_mem[UP_G] * b_sqr + gas_idx * n_phases] += pp[i + gas_idx * n_phases];
00690 if (is_o) reg_values_ [i + m_mem[UP_O] * b_sqr + oil_idx * n_phases] += pp[i + oil_idx * n_phases];
00691 if (is_w) reg_values_ [i + m_mem[UP_W] * b_sqr + wat_idx * n_phases] += pp[i + wat_idx * n_phases];
00692
00693 if (is_g) reg_values_ [i + p_mem[UP_G] * b_sqr + gas_idx * n_phases] -= pp[i + gas_idx * n_phases];
00694 if (is_o) reg_values_ [i + p_mem[UP_O] * b_sqr + oil_idx * n_phases] -= pp[i + oil_idx * n_phases];
00695 if (is_w) reg_values_ [i + p_mem[UP_W] * b_sqr + wat_idx * n_phases] -= pp[i + wat_idx * n_phases];
00696 }
00697 }
00698
00707 struct cfl_info
00708 {
00709 rhs_item_array_t f11_i;
00710 rhs_item_array_t f12_i;
00711 rhs_item_array_t f21_i;
00712 rhs_item_array_t f22_i;
00713
00714 cfl_info (index_t n_elems)
00715 {
00716 f11_i.assign (n_elems, 0.0f);
00717 f12_i.assign (n_elems, 0.0f);
00718 f21_i.assign (n_elems, 0.0f);
00719 f22_i.assign (n_elems, 0.0f);
00720 }
00721
00722 void clear (index_t n_elems)
00723 {
00724 f11_i.assign (n_elems, 0.0f);
00725 f12_i.assign (n_elems, 0.0f);
00726 f21_i.assign (n_elems, 0.0f);
00727 f22_i.assign (n_elems, 0.0f);
00728 }
00729 rhs_item_t calc_F_i (item_t dt, item_t i)
00730 {
00731 return dt*fabs(0.5*(f11_i[i]+f22_i[i]+sqrt((f11_i[i]+f22_i[i])*(f11_i[i]+f22_i[i])-4*(f11_i[i]*f22_i[i]-f12_i[i]*f21_i[i]))));
00732 }
00733 };
00734
00735
00745 void
00746 fill_cfl (cfl_info &f, const rhs_item_t *truns, index_t i_cell, index_t j_cell, rhs_item_array_t &cfl, item_array_t &saturation_3p)
00747 {
00748
00749 int k;
00750 index_t i, j, k_cell;
00751 const data_t &data_ic = data_[i_cell];
00752 const data_t &data_jc = data_[j_cell];
00753
00754 item_t mobility_tot = 0;
00755 if (is_g) mobility_tot += MOBILITY_UP_G;
00756 if (is_o) mobility_tot += MOBILITY_UP_O;
00757 if (is_w) mobility_tot += MOBILITY_UP_W;
00758
00759 item_t p_cwoi_deriv = data_[i_cell].s_deriv_cap_pressure[ds_w];
00760 item_t p_cwoj_deriv = data_[j_cell].s_deriv_cap_pressure[ds_w];
00761 item_t p_cgoi_deriv = 0;
00762 item_t p_cgoj_deriv = 0;
00763 item_t p_cwoi = data_[i_cell].cap_pressure[ds_w];
00764 item_t p_cwoj = data_[j_cell].cap_pressure[ds_w];
00765 item_t p_cgoi = 0;
00766 item_t p_cgoj = 0;
00767
00768
00769 if (is_g == 1)
00770 {
00771 p_cgoi_deriv = data_[i_cell].s_deriv_cap_pressure[ds_g];
00772 p_cgoj_deriv = data_[j_cell].s_deriv_cap_pressure[ds_g];
00773 p_cgoi = data_[i_cell].cap_pressure[ds_g];
00774 p_cgoj = data_[j_cell].cap_pressure[ds_g];
00775 }
00776
00777
00778 item_t delta_psi_w = fabs(psi_[d_w]);
00779 item_t delta_psi_o = fabs(psi_[d_o]);
00780
00781 item_t delta_psi_g;
00782 if (is_g == 1) delta_psi_g = fabs(psi_[d_g]);
00783
00784
00785 if (is_g == 0)
00786 {
00787 f.f11_i[i_cell] += (MOBILITY_UP_O*SW_DERIV_MOBILITY_UP_W*delta_psi_w
00788 - MOBILITY_UP_W*SW_DERIV_MOBILITY_UP_O*delta_psi_o - truns[0]*MOBILITY_UP_W*MOBILITY_UP_O
00789 *(p_cwoi_deriv + p_cwoj_deriv))/mobility_tot;
00790
00791 f.f11_i[j_cell] = f.f11_i[i_cell];
00792 }
00793 else
00794 {
00795 f.f11_i[i_cell] += ((MOBILITY_UP_G + MOBILITY_UP_O)*SW_DERIV_MOBILITY_UP_W*delta_psi_w
00796 - MOBILITY_UP_W*SW_DERIV_MOBILITY_UP_O*delta_psi_o - truns[0]*MOBILITY_UP_W*(MOBILITY_UP_G + MOBILITY_UP_O)
00797 *(p_cwoi_deriv + p_cwoj_deriv))/mobility_tot;
00798
00799 f.f11_i[j_cell] = f.f11_i[i_cell];
00800 }
00801
00802 if (is_g == 0)
00803 {
00804 f.f12_i[i_cell] = 0;
00805 f.f12_i[j_cell] = f.f12_i[i_cell];
00806 }
00807
00808 else
00809 {
00810 f.f12_i[i_cell] += -((MOBILITY_UP_W*SG_DERIV_MOBILITY_UP_O*delta_psi_o) + MOBILITY_UP_W*SG_DERIV_MOBILITY_UP_G
00811 *delta_psi_g-(MOBILITY_UP_O + MOBILITY_UP_G)*SG_DERIV_MOBILITY_UP_W*delta_psi_w + truns[0]*MOBILITY_UP_W*MOBILITY_UP_G
00812 *(p_cgoi_deriv + p_cgoj_deriv))/mobility_tot;
00813
00814 f.f12_i[j_cell] = f.f12_i[i_cell];
00815 }
00816
00817 if (is_g == 0)
00818 {
00819 f.f21_i[i_cell] = 0;
00820 f.f21_i[j_cell] = f.f21_i[i_cell];
00821 }
00822
00823 else
00824 {
00825 f.f21_i[i_cell] += -((MOBILITY_UP_G*SW_DERIV_MOBILITY_UP_W)*delta_psi_w + MOBILITY_UP_G*SW_DERIV_MOBILITY_UP_O
00826 *delta_psi_o - truns[0]*MOBILITY_UP_G*MOBILITY_UP_W*(p_cgoi_deriv + p_cgoj_deriv))/mobility_tot;
00827
00828 f.f21_i[j_cell] = f.f21_i[i_cell];
00829
00830 }
00831
00832 if (is_g == 0)
00833 {
00834 f.f22_i[i_cell] = 0;
00835 f.f22_i[j_cell] = f.f22_i[i_cell];
00836 }
00837
00838 else
00839 {
00840 f.f22_i[i_cell] += (-MOBILITY_UP_G*SG_DERIV_MOBILITY_UP_O*delta_psi_o + (MOBILITY_UP_W + MOBILITY_UP_O)
00841 *SG_DERIV_MOBILITY_UP_G*delta_psi_g - MOBILITY_UP_G*SG_DERIV_MOBILITY_UP_W*delta_psi_w
00842 + truns[0]*MOBILITY_UP_G*(MOBILITY_UP_O + MOBILITY_UP_W)*(p_cgoi_deriv + p_cgoj_deriv))/mobility_tot;
00843
00844 f.f22_i[j_cell] = f.f22_i[i_cell];
00845 }
00846
00847
00848 if (fabs (saturation_3p[i_cell * n_phases + d_w] - saturation_3p[j_cell * n_phases + d_w]) > 0.0001 ||
00849 fabs (saturation_3p[i_cell * n_phases + d_o] - saturation_3p[j_cell * n_phases + d_o]) > 0.0001)
00850 {
00851 cfl[i_cell] = UNKNOWN;
00852 cfl[j_cell] = UNKNOWN;
00853 }
00854 }
00855
00856 public:
00857
00858 boost::array <item_t, b_sqr + n_phases> pp;
00859
00860 index_t d_w;
00861 index_t d_g;
00862 index_t d_o;
00863 index_t ds_w;
00864 index_t ds_g;
00865 const item_array_t &saturation_3p_;
00866 const item_array_t &pressure_;
00867 const item_array_t &gas_oil_ratio_;
00868 const main_var_array_t &main_vars_;
00869 const data_array_t &data_;
00870 const item_array_t &depths_;
00871 item_t gravity_;
00872 index_t n_sec_vars;
00873
00874 const rhs_item_array_t &sp_diag_;
00875 const rhs_item_array_t &s_rhs_;
00876 rhs_item_array_t ®_values_;
00877 rhs_item_array_t &rhs_;
00878
00879 boost::array <item_t, n_phases> rho_;
00880 boost::array <item_t, n_phases> psi_;
00881 boost::array <index_t, n_phases> up_;
00882 boost::array <index_t, n_phases> up_cell_;
00883 boost::array <item_t, n_phases> flow_;
00884
00885 const index_t *m_mem;
00886 const index_t *p_mem;
00887
00888 item_t depth_top;
00889 item_t sw_deriv_mobility_up_w;
00890 item_t sw_deriv_mobility_up_g;
00891 item_t sw_deriv_mobility_up_o;
00892 item_t sg_deriv_mobility_up_w;
00893 item_t sg_deriv_mobility_up_g;
00894 item_t sg_deriv_mobility_up_o;
00895 item_t so_deriv_mobility_up_w;
00896 item_t so_deriv_mobility_up_g;
00897 item_t so_deriv_mobility_up_o;
00898 item_t p_deriv_mobility_up_w;
00899 item_t p_deriv_mobility_up_g;
00900 item_t p_deriv_mobility_up_o;
00901 item_t mobility_up_w;
00902 item_t mobility_up_g;
00903 item_t mobility_up_o;
00904 item_t gor_up_o;
00905 item_t p_deriv_gor_up_o;
00906 main_var_type main_var_up_g;
00907 main_var_type main_var_up_o;
00908 };
00909
00910 template <typename fi_operator_t, typename mpfa_t>
00911 void
00912 block_connections_mpfa (typename mpfa_t::item_t dt, fi_operator_t &fi_operator, mpfa_t mpfa_)
00913 {
00914 typedef typename mpfa_t::index_t index_t;
00915 typedef typename mpfa_t::item_t item_t;
00916 typedef typename mpfa_t::rhs_item_t rhs_item_t;
00917 typedef typename mpfa_t::item_array_t item_array_t;
00918
00919 #ifdef BS_BOS_CORE_USE_CSR_ILU_CFL_PREC
00920 index_t n_active_cells = fi_operator.n_cells_;
00921 bool use_cfl = fi_operator.calc_model_->ts_params->get_bool (fi_params::USE_CFL);
00922 static typename mpfa_t::cfl_info cfl (n_active_cells);
00923
00924 if (use_cfl)
00925 {
00926 cfl.clear (n_active_cells);
00927 assign (fi_operator.cfl_, n_active_cells, 0);
00928 }
00929 #endif
00930
00931
00932
00933 for (index_t con = 0; con < fi_operator.n_connections_; ++con)
00934 {
00935 const index_t row[] = {fi_operator.trns_rows_ptr_[con], fi_operator.trns_rows_ptr_[con + 1]};
00936 BS_ASSERT (row[1] - row[0] >= 2) (row[0]) (row[1]);
00937
00938 const index_t *cells = &fi_operator.trns_cols_ptr_[row[0]];
00939 const rhs_item_t *truns = &fi_operator.trns_values_[row[0]];
00940
00941 index_t i = cells[0];
00942 index_t j = cells[1];
00943
00944 mpfa_.m_mem = &fi_operator.m_array_[row[0]];
00945 mpfa_.p_mem = &fi_operator.p_array_[row[0]];
00946
00947 mpfa_.compute_rho (cells, row[1] - row[0]);
00948 mpfa_.compute_psi (truns, cells, row[1] - row[0]);
00949 mpfa_.compute_upstream (0, 1, i, j);
00950 mpfa_.compute_flow (dt);
00951
00952 mpfa_.fill_rhs (i, j);
00953
00954 for (index_t k = 0, kcnt = row[1] - row[0]; k < kcnt; ++k)
00955 {
00956 mpfa_.compute_k_derivs (k, cells[k], fi_operator.main_vars_[cells[k]], dt * truns [k], i, j);
00957 }
00958
00959 mpfa_.compute_up_derivs (dt, i, j);
00960
00961
00962 #ifdef BS_BOS_CORE_USE_CSR_ILU_CFL_PREC
00963 if (use_cfl)
00964 {
00965
00966 mpfa_.fill_cfl (cfl, truns, i, j, fi_operator.cfl_, fi_operator.saturation_3p_);
00967 }
00968 #endif
00969 }
00970
00971 #ifdef BS_BOS_CORE_USE_CSR_ILU_CFL_PREC
00972 if (use_cfl)
00973 {
00974 for (index_t k = 0; k < n_active_cells; ++k)
00975 {
00976 if (fi_operator.cfl_[k] == UNKNOWN)
00977 {
00978 fi_operator.cfl_[k] = cfl.calc_F_i (dt, k) / fi_operator.volume_[k];
00979 }
00980 else
00981 fi_operator.cfl_[k] = IS_STABLE;
00982 }
00983 }
00984 #endif
00985 }
00986 }
00987
00993 template <typename strategy_t, bool is_w, bool is_g, bool is_o>
00994 bool
00995 fi_operator_impl <strategy_t, is_w, is_g, is_o>::block_connections_mpfa (const item_t &dt)
00996 {
00997 mpfa::block_connections_mpfa (dt, *this, mpfa::mpfa_base_impl <strategy_t, is_w, is_g, is_o> (*this, flux_rhs_, reg_values_));
00998 return false;
00999 }
01000
01001
01002
01003 }
01004
01005
01006
01007 #endif // BS_FI_OPERATOR_BLOCK_CONNECTION_MPFA_2_H_
01008
