TheBSEagleProject: fi_operator_block_connections

00001 
00010 #ifndef BS_FI_OPERATOR_BLOCK_CONNECTIONS_MPFA_H_
00011 #define BS_FI_OPERATOR_BLOCK_CONNECTIONS_MPFA_H_
00012 
00013 namespace blue_sky {
00014 
00015   namespace tpfa 
00016   {
00017     template <typename strategy_t>
00018     struct mpfa_impl 
00019     {
00020     public:
00021       typedef typename strategy_t::item_t             item_t;
00022       typedef typename strategy_t::index_t            index_t;
00023       typedef typename strategy_t::item_array_t       item_array_t;
00024       typedef typename strategy_t::index_array_t      index_array_t;
00025       typedef boost::array <index_t, FI_PHASE_TOT>    up_cell_array_t;
00026       typedef calc_model <strategy_t>                 calc_model_t;
00027       typedef typename calc_model_t::data_t           data_t;
00028       typedef typename calc_model_t::data_array_t     data_array_t;
00029       typedef typename calc_model_t::main_var_array_t main_var_array_t;
00030 
00031     public:
00032       mpfa_impl (const fi_operator_impl <strategy_t> &fi_operator_, item_array_t &rhs_, item_array_t &reg_values_)
00033       : fi_operator_ (fi_operator_),
00034       n_phases (fi_operator_.n_phases),
00035       is_w (fi_operator_.is_w),
00036       is_g (fi_operator_.is_g),
00037       is_o (fi_operator_.is_o),
00038       d_w (fi_operator_.d_w),
00039       d_g (fi_operator_.d_g),
00040       d_o (fi_operator_.d_o),
00041       saturation_3p_ (fi_operator_.saturation_3p_),
00042       pressure_ (fi_operator_.pressure_),
00043       gas_oil_ratio_ (fi_operator_.gas_oil_ratio_),
00044       main_vars_ (fi_operator_.main_vars_),
00045       data_ (fi_operator_.data_),
00046       gravity_ (fi_operator_.calc_model_->internal_constants.gravity_constant),
00047       n_sec_vars (fi_operator_.n_sec_vars),
00048       rhs_ (rhs_),
00049       reg_values_ (reg_values_)
00050       {
00051       }
00052 
00053       /*BS_FORCE_INLINE*/ void
00054       mpfa_calc_avarage_density(item_t *ave_density,
00055           const index_t *cell_ind_block,
00056           const int &n_cells_in_conn) const
00057       {
00058         item_t counter[FI_PHASE_TOT] = {0};
00059         index_t flag[FI_PHASE_TOT] = {0};
00060 
00061         memset (ave_density, 0, sizeof (item_t) * FI_PHASE_TOT);
00062         //ave_density[0] = ave_density[1] = ave_density[2] = 0;
00063 
00064         for (int j = 0; j < n_cells_in_conn; ++j)
00065           {
00066             index_t ic = cell_ind_block[j];
00067             index_t idx = ic * n_phases;
00068             const data_t &data_ic = data_[ic];
00069 
00070             if (n_phases > 1)
00071               {
00072                 bool is_w_ = is_w && saturation_3p_[idx + d_w] > EPS_DIV;
00073                 bool is_g_ = is_g && saturation_3p_[idx + d_g] > EPS_DIV;
00074                 bool is_o_ = is_o && saturation_3p_[idx + d_o] > EPS_DIV;
00075 
00076                 if (is_w_)
00077                   {
00078                     ave_density[d_w] += data_ic.density[d_w];
00079                     ++counter[d_w];
00080                     ++flag[d_w];
00081                   }
00082                 if (is_g_)
00083                   {
00084                     ave_density[d_g] += data_ic.density[d_g];
00085                     ++counter[d_g];
00086                     ++flag[d_g];
00087                   }
00088                 if (is_o_)
00089                   {
00090                     ave_density[d_o] += data_ic.density[d_o];
00091                     ++counter[d_o];
00092                     ++flag[d_o];
00093                   }
00094               }
00095             else
00096               {
00097                 ave_density[0] += data_ic.density[0];
00098                 ++counter[0];
00099                 ++flag[0];
00100               }
00101           }
00102         for (int j = 0; j < n_phases; ++j)
00103           {
00104             if (flag[j])
00105               ave_density[j] /= counter[j];
00106           }
00107       }
00108 
00109       BS_FORCE_INLINE void
00110       mpfa_calc_potential (item_t *cell_pot,
00111           item_t *sum_cell_pot,
00112           item_t *ave_density,
00113           const item_t *truns_block,
00114           const item_array_t &depth,
00115           const index_t * cell_ind_block,
00116           index_t cell_m,
00117           item_t dt,
00118           int n_cells_in_conn) const
00119       {
00120         item_t g_h      = gravity_;
00121         item_t base_h   = depth[cell_m]; 
00122         sum_cell_pot[0] = sum_cell_pot[1] = sum_cell_pot[2] = 0;
00123 
00124         for (index_t j = 0; j < n_cells_in_conn; ++j)
00125           {
00126             index_t ic            = cell_ind_block[j];
00127             item_t dh             = depth[ic] - base_h; 
00128             const data_t &data_ic = data_[ic];
00129             item_t truns_block_dt = truns_block[j] * dt;
00130 
00131             if (is_w)
00132               {
00133                 index_t k = j * n_phases + d_w;
00134                 cell_pot[k] = pressure_[ic] - ave_density[d_w] * g_h * dh;
00135                 if (this->n_phases > 1)
00136                   cell_pot[k] += data_ic.cap_pressure[d_w];
00137                 sum_cell_pot[d_w] += cell_pot[k] * truns_block_dt;
00138               }
00139             if (is_g)
00140               {
00141                 index_t k = j * this->n_phases + d_g;
00142                 cell_pot[k] = pressure_[ic] - ave_density[d_g] * g_h * dh;
00143                 if (this->n_phases > 1)
00144                   cell_pot[k] += data_ic.cap_pressure[d_g];
00145                 sum_cell_pot[d_g] += cell_pot[k] * truns_block_dt;
00146               }
00147             if (is_o)
00148               {
00149                 index_t k = j * this->n_phases + d_o;
00150                 cell_pot[k] = pressure_[ic] - ave_density[d_o] * g_h * dh;
00151                 sum_cell_pot[d_o] += cell_pot[k] * truns_block_dt;
00152               }
00153           }
00154       }
00155 
00156       /*BS_FORCE_INLINE */void
00157       mpfa_calc_upstream(up_cell_array_t &up_cell,
00158           item_t *sum_cell_pot,
00159           const index_t &cell_m,
00160           const index_t &cell_p) const
00161       {
00162         index_t cell_[] = {cell_m, cell_p};
00163         for (int j = 0; j < n_phases; ++j)
00164           {
00165             up_cell[j] = cell_[sum_cell_pot[j] > 0 ? 0 : 1];
00166           }
00167       }
00168 
00169       /*BS_FORCE_INLINE */void
00170       mpfa_fill_rhs(item_t * & rhs_m,
00171                                             item_t * & rhs_p,
00172                                             const up_cell_array_t &up_cell,
00173                                             const item_t *truns_block,
00174                                             item_t *cell_pot,
00175                                             item_array_t &flux_rhs,
00176                                             const item_t &dt,
00177                                             const int &n_cells_in_conn,
00178                                             const index_t &cell_m,
00179                                             const index_t &cell_p,
00180                                             const index_t &equ_w, const index_t equ_g, const index_t &equ_o) const
00181       {
00182         rhs_m = &flux_rhs[0] + cell_m * this->n_phases;
00183         rhs_p = &flux_rhs[0] + cell_p * this->n_phases;
00184         item_t rhs_block_m[3] = {0}, rhs_block_p[3] = {0};
00185         //rhs_block_m[0] = rhs_block_m[1] = rhs_block_m[2] = 0;
00186         //rhs_block_p[0] = rhs_block_p[1] = rhs_block_p[2] = 0;
00187 
00188         static data_t dummy_data;
00189         const data_t &data_m = data_[cell_m];
00190         const data_t &data_w = is_w ? data_[up_cell[d_w]] : dummy_data;
00191         const data_t &data_g = is_g ? data_[up_cell[d_g]] : dummy_data;
00192         const data_t &data_o = is_o ? data_[up_cell[d_o]] : dummy_data;
00193 
00194         item_t flow_m = data_m.truns_mult * dt;
00195         item_t flow_w = is_w ? flow_m * data_w.mobility[d_w] : 0;
00196         item_t flow_g = is_g ? flow_m * data_g.mobility[d_g] : 0;
00197         item_t flow_o = is_o ? flow_m * data_o.mobility[d_o] : 0;
00198 
00199         index_t k_j = 0;
00200         for (int j = 0; j < n_cells_in_conn; ++j, k_j += n_phases)
00201           {
00202             // water
00203             if (is_w)
00204               {
00205                 index_t k = k_j + d_w;
00206                 item_t flow = flow_w * truns_block[j] * cell_pot[k];
00207                 rhs_block_m[equ_w] -= flow;
00208                 rhs_block_p[equ_w] += flow;
00209               }
00210             // oil
00211             if (is_o)
00212               {
00213                 index_t k = k_j + d_o;
00214                 item_t flow = flow_o * truns_block[j] * cell_pot[k];
00215                 rhs_block_m[equ_o] -= flow;
00216                 rhs_block_p[equ_o] += flow;
00217                 // disolved gas
00218                 if (is_g)
00219                   {
00220                     flow *= gas_oil_ratio_[up_cell[d_o]];
00221                     rhs_block_m[equ_g] -= flow;
00222                     rhs_block_p[equ_g] += flow;
00223                   }
00224               }
00225             // gas
00226             if (is_g)
00227               {
00228                 index_t k = k_j + d_g;
00229                 item_t flow = flow_g * truns_block[j] * cell_pot[k];
00230                 rhs_block_m[equ_g] -= flow;
00231                 rhs_block_p[equ_g] += flow;
00232               }
00233           }
00234         for (int j = 0; j < n_phases; ++j)
00235           {
00236             rhs_m[j] += rhs_block_m[j];
00237             rhs_p[j] += rhs_block_p[j];
00238           }
00239       }
00240 
00241       /*BS_FORCE_INLINE*/ void
00242       mpfa_fill_jacobian( boost::array <item_t, 18>  &m_mem,
00243           boost::array <item_t, 18>  &p_mem,
00244           const up_cell_array_t &up_cell,
00245           const item_t *cell_pot,
00246           const item_t *truns_block,
00247           item_t * &rhs_m,
00248           item_t * &rhs_p,
00249           const index_t *cell_ind_block,
00250           item_array_t &sp_diag,
00251           item_array_t &s_rhs,
00252           const item_t &dt,
00253           const int &n_cells_in_conn,
00254           const index_t &cell_m,
00255           const index_t &cell_p,
00256           const index_t &equ_w, const index_t equ_g, const index_t &equ_o) const
00257       {
00258         //m_mem.assign(0);
00259         //p_mem.assign(0);
00260         assign (m_mem, item_t (0));
00261         assign (p_mem, item_t (0));
00262 
00263         static data_t dummy_data;
00264 
00265         item_t *sp_diag_val = &sp_diag[0];
00266         item_t *s_rhs_val = &s_rhs[0];
00267 
00268         const data_t *data_m = &data_[cell_m];
00269         const data_t *data_w = &dummy_data;
00270         const data_t *data_o = &dummy_data;
00271         const data_t *data_g = &dummy_data;
00272 
00273         index_t up_cell_dw  = 0;
00274         index_t up_cell_do  = 0;
00275         index_t up_cell_dg  = 0;
00276 
00277         index_t sp_size_w   = 0;
00278         index_t sp_size_o   = 0;
00279         index_t sp_size_g   = 0;
00280 
00281         index_t rhs_size_w  = 0;
00282         index_t rhs_size_o  = 0;
00283         index_t rhs_size_g  = 0;
00284 
00285         item_t truns_mob_w  = 0;
00286         item_t truns_mob_o  = 0;
00287         item_t truns_mob_g  = 0;
00288 
00289         item_t truns_p_mob_w  = 0;
00290         item_t truns_p_mob_o  = 0;
00291         item_t truns_p_mob_g  = 0;
00292 
00293         item_t truns_s_mob_ww  = 0;
00294         item_t truns_s_mob_ow  = 0;
00295         item_t truns_s_mob_og  = 0;
00296         item_t truns_s_mob_oo  = 0;
00297         item_t truns_s_mob_gg  = 0;
00298 
00299         // calculate mobility indexies
00300         index_t mob_w_dw = d_w * n_phases + d_w;
00301         index_t mob_o_dw = d_o * n_phases + d_w;
00302         index_t mob_o_dg = d_o * n_phases + d_g;
00303         index_t mob_o_do = d_o * n_phases + d_o;
00304         index_t mob_g_dg = d_g * n_phases + d_g;
00305 
00306         if (is_w)
00307           {
00308             data_w      = &data_[up_cell[d_w]];
00309             up_cell_dw  = up_cell[d_w];
00310             sp_size_w   = n_sec_vars * this->n_phases * up_cell_dw;
00311             rhs_size_w  = n_sec_vars * up_cell_dw;
00312             truns_mob_w = data_m->truns_mult * dt * data_w->mobility[d_w];
00313             truns_p_mob_w = data_m->truns_mult *dt * data_w->p_deriv_mobility[d_w];
00314             truns_s_mob_ww  = data_m->truns_mult * dt * data_w->s_deriv_mobility[mob_w_dw];
00315          }
00316 
00317         if (is_o)
00318           {
00319             data_o      = &data_[up_cell[d_o]];
00320             up_cell_do  = up_cell[d_o];
00321             rhs_size_o  = n_sec_vars * up_cell_do;
00322             sp_size_o   = n_sec_vars * this->n_phases * up_cell_do;
00323             truns_mob_o = data_m->truns_mult * dt * data_o->mobility[d_o];
00324             truns_p_mob_o = data_m->truns_mult *dt * data_o->p_deriv_mobility[d_o];
00325             truns_s_mob_ow  = data_m->truns_mult * dt * data_o->s_deriv_mobility[mob_o_dw];
00326             truns_s_mob_og  = data_m->truns_mult * dt * data_o->s_deriv_mobility[mob_o_dg];
00327             truns_s_mob_oo  = data_m->truns_mult * dt * data_o->s_deriv_mobility[mob_o_do];
00328           }
00329         if (is_g)
00330           {
00331             data_g      = &data_[up_cell[d_g]];
00332             up_cell_dg  = up_cell[d_g];
00333             sp_size_g   = n_sec_vars * this->n_phases * up_cell_dg;
00334             rhs_size_g  = n_sec_vars * up_cell_dg;
00335             truns_mob_g = data_m->truns_mult * dt * data_g->mobility[d_g];
00336             truns_p_mob_g = data_m->truns_mult *dt * data_g->p_deriv_mobility[d_g];
00337             truns_s_mob_gg  = data_m->truns_mult * data_g->s_deriv_mobility[mob_g_dg];
00338          }
00339 
00340         index_t k_j = 0;
00341         for (int j = 0; j < n_cells_in_conn; ++j, k_j += n_phases)
00342           {
00343             index_t ic = cell_ind_block[j];
00344 
00345             index_t k;
00346             item_t dsw, dsg, dso, dpo;
00347             item_t dsw_up, dsg_up, dpo_up, dso_up;
00348 
00349             item_t *pp_block_m = &m_mem [j * 9];
00350             item_t *pp_block_p = &p_mem [j * 9];
00351             item_t *pp_block_m_up;
00352             item_t *pp_block_p_up;
00353             item_t ps_block_p_up[3];
00354             item_t ps_block_m_up[3];
00355             item_t *sp_block_up;
00356             item_t *s_rhs_block_up;
00357 
00358             item_t ps_block_m[FI_PHASE_TOT] = {0};   // temporary block for PS
00359             item_t ps_block_p[FI_PHASE_TOT] = {0};   // temporary block for PS
00360 
00361             item_t *sp_block = sp_diag_val + n_sec_vars * this->n_phases * ic;
00362             item_t *s_rhs_block = s_rhs_val + n_sec_vars * ic;
00363 
00364             item_t truns_block_j = truns_block[j];
00365 
00366             // water
00367             if (is_w)
00368               {
00369                 ps_block_p_up[0]  = ps_block_p_up[1] = ps_block_p_up[2] = 0;
00370                 ps_block_m_up[0]  = ps_block_m_up[1] = ps_block_m_up[2] = 0;
00371                 sp_block_up       = sp_diag_val + sp_size_w;
00372                 s_rhs_block_up    = s_rhs_val + rhs_size_w;
00373 
00374                 if (up_cell_dw == cell_ind_block[0])
00375                   {
00376                     pp_block_m_up = &m_mem[0];
00377                     pp_block_p_up = &p_mem[0];
00378                   }
00379                 else
00380                   {
00381                     pp_block_m_up = &m_mem[9];
00382                     pp_block_p_up = &p_mem[9];
00383                   }
00384                 k = k_j + d_w;
00385                 item_t cell_pot_k = cell_pot[k];
00386                 dpo = truns_mob_w * truns_block[j];
00387                 dpo_up = 0;
00388                 if (ic == up_cell_dw)
00389                   {
00390                     dpo += truns_p_mob_w * truns_block_j * cell_pot_k;
00391                   }
00392                 else
00393                   {
00394                     dpo_up = truns_p_mob_w * truns_block_j * cell_pot_k;
00395                   }
00396 
00397                 if (n_phases > 1)
00398                   {
00399                     dsg = 0;
00400                     dso = 0;
00401                     dsw = truns_mob_w * truns_block_j * data_[ic].s_deriv_cap_pressure[d_w];
00402                     dsg_up = dso_up = dsw_up = 0;
00403                     if (ic == up_cell_dw)
00404                       {
00405                         dsw += truns_s_mob_ww * truns_block_j * cell_pot_k;
00406                       }
00407                     else
00408                       {
00409                         dsw_up = truns_s_mob_ww * truns_block_j * cell_pot_k;
00410                       }
00411                     if (n_phases == 3)
00412                       {
00413                         pp_block_m[p3_wat_po]     += dpo;
00414                         pp_block_m_up[p3_wat_po]  += dpo_up;
00415                         pp_block_p[p3_wat_po]     -= dpo;
00416                         pp_block_p_up[p3_wat_po]  -= dpo_up;
00417 
00418                         ps_block_m[p3_wat]        += dsw;
00419                         ps_block_m_up[p3_wat]     += dsw_up;
00420                         ps_block_p[p3_wat]        += -dsw;
00421                         ps_block_p_up[p3_wat]     += -dsw_up;
00422                       }
00423                     else
00424                       {
00425                         pp_block_m[p2ow_wat_po]     += dpo;
00426                         pp_block_m_up[p2ow_wat_po]  += dpo_up;
00427                         pp_block_p[p2ow_wat_po]     -= dpo;
00428                         pp_block_p_up[p2ow_wat_po]  -= dpo_up;
00429 
00430                         ps_block_m[p2ow_wat]        += dsw;
00431                         ps_block_m_up[p2ow_wat]     += dsw_up;
00432                         ps_block_p[p2ow_wat]        += -dsw;
00433                         ps_block_p_up[p2ow_wat]     += -dsw_up;
00434                       }
00435 
00436                     M_MINUS_VV_PROD (n_phases, ps_block_m_up, sp_block_up, pp_block_m_up);
00437                     V_MINUS_VS_PROD (n_phases, ps_block_m_up, s_rhs_block_up, rhs_m);
00438 
00439                     M_MINUS_VV_PROD (n_phases, ps_block_p_up, sp_block_up, pp_block_p_up);
00440                     V_MINUS_VS_PROD (n_phases, ps_block_p_up, s_rhs_block_up, rhs_p);
00441                   }
00442                 else
00443                   {
00444                     pp_block_m[0]     += dpo;
00445                     pp_block_m_up[0]  += dpo_up;
00446                     pp_block_p[0]     -= dpo;
00447                     pp_block_p_up[0]  -= dpo_up;
00448                   }
00449               }
00450             // oil
00451             if (is_o)
00452               {
00453                 item_t g_dpo, g_dso, g_dsw, g_dsg;
00454                 item_t g_dpo_up, g_dso_up, g_dsw_up, g_dsg_up;
00455                 k = k_j + d_o;
00456                 item_t cell_pot_k = cell_pot[k];
00457 
00458                 ps_block_p_up[0] = ps_block_p_up[1] = ps_block_p_up[2] = 0;
00459                 ps_block_m_up[0] = ps_block_m_up[1] = ps_block_m_up[2] = 0;
00460                 sp_block_up = sp_diag_val + sp_size_o;
00461                 s_rhs_block_up = s_rhs_val + rhs_size_o;
00462 
00463                 if (up_cell_do == cell_ind_block[0])
00464                   {
00465                     pp_block_m_up = &m_mem[0];
00466                     pp_block_p_up = &p_mem[0];
00467                   }
00468                 else
00469                   {
00470                     pp_block_m_up = &m_mem[9];
00471                     pp_block_p_up = &p_mem[9];
00472                   }
00473 
00474                 // calculate pressure derivatives
00475                 dpo = truns_mob_o * truns_block_j;
00476                 dpo_up = 0;
00477                 if (ic == up_cell_do)
00478                   {
00479                     dpo += truns_p_mob_o * truns_block_j * cell_pot_k;
00480                   }
00481                 else
00482                   {
00483                     dpo_up = truns_p_mob_o * truns_block_j * cell_pot_k;
00484                   }
00485                 if (is_g)
00486                   {
00487                     g_dpo = gas_oil_ratio_[up_cell_do] * dpo;
00488                     g_dpo_up = gas_oil_ratio_[up_cell_do] * dpo_up;
00489                     if (ic == up_cell_do && FI_CHK_SG (main_vars_, ic))
00490                       {
00491                         g_dpo += truns_mob_o * truns_block_j * data_o->p_deriv_gas_oil_ratio * cell_pot_k;
00492                       }
00493                     else if (FI_CHK_SG (main_vars_, up_cell_do))
00494                       {
00495                         g_dpo_up += truns_mob_o * truns_block_j * data_o->p_deriv_gas_oil_ratio * cell_pot_k;
00496                       }
00497                   }
00498 
00499                 dsw = 0;
00500                 dso = 0;
00501                 dsg = 0;
00502                 dsw_up = 0;
00503                 dso_up = 0;
00504                 dsg_up = 0;
00505 
00506                 // calculate saturation deriv
00507                 if (ic == up_cell[d_o])
00508                   {
00509                     if (is_w)
00510                       dsw = truns_s_mob_ow * truns_block_j * cell_pot_k;
00511                     if (is_g)
00512                       dsg = truns_s_mob_og * truns_block_j * cell_pot_k;
00513                     if (n_phases > 1)
00514                       dso = truns_s_mob_oo * truns_block_j * cell_pot_k;
00515                   }
00516                 else
00517                   {
00518                     if (is_w)
00519                       dsw_up = truns_s_mob_ow * truns_block_j * cell_pot_k;
00520                     if (is_g)
00521                       dsg_up = truns_s_mob_og * truns_block_j * cell_pot_k;
00522                     if (n_phases > 1)
00523                       dso_up = truns_s_mob_oo * truns_block_j * cell_pot_k;
00524                   }
00525                 if (is_g)
00526                   {
00527                     g_dso = gas_oil_ratio_[up_cell_do] * dso;
00528                     g_dso_up = gas_oil_ratio_[up_cell_do] * dso_up;
00529                     if (is_w)
00530                       {
00531                         g_dsw = gas_oil_ratio_[up_cell_do] * dsw;
00532                         g_dsw_up = gas_oil_ratio_[up_cell_do] * dsw_up;
00533                       }
00534 
00535                     g_dsg = gas_oil_ratio_[up_cell_do] * dsg;
00536                     g_dsg_up = gas_oil_ratio_[up_cell_do] * dsg_up;
00537 
00538                     if (ic == up_cell[d_o] && FI_CHK_RO (main_vars_, ic))
00539                       {
00540                         g_dsg += truns_s_mob_og * truns_block_j * cell_pot_k;
00541                       }
00542                     else if (FI_CHK_RO (main_vars_, ic))
00543                       {
00544                         g_dsg_up += truns_s_mob_og * truns_block_j * cell_pot_k;
00545                       }
00546                   }
00547 
00548                 if (n_phases == 3)
00549                   {
00550                     //oil
00551                     pp_block_m[p3_oil_po] += dpo;
00552                     pp_block_m[p3_oil_so] += dso;
00553                     pp_block_m[p3_oil_sg] += dsg;
00554                     ps_block_m[p3_oil] += dsw;
00555 
00556                     pp_block_p[p3_oil_po] -= dpo;
00557                     pp_block_p[p3_oil_so] -= dso;
00558                     pp_block_p[p3_oil_sg] -= dsg;
00559                     ps_block_p[p3_oil] += -dsw;
00560 
00561                     pp_block_m_up[p3_oil_po] += dpo_up;
00562                     pp_block_m_up[p3_oil_so] += dso_up;
00563                     pp_block_m_up[p3_oil_sg] += dsg_up;
00564                     ps_block_m_up[p3_oil] += dsw_up;
00565 
00566                     pp_block_p_up[p3_oil_po] -= dpo_up;
00567                     pp_block_p_up[p3_oil_so] -= dso_up;
00568                     pp_block_p_up[p3_oil_sg] -= dsg_up;
00569                     ps_block_p_up[p3_oil] += -dsw_up;
00570 
00571                     // gas
00572                     pp_block_m[p3_gas_po] += g_dpo;
00573                     pp_block_m[p3_gas_so] += g_dso;
00574                     pp_block_m[p3_gas_sg] += g_dsg;
00575                     ps_block_m[p3_gas] += g_dsw;
00576 
00577                     pp_block_p[p3_gas_po] -= g_dpo;
00578                     pp_block_p[p3_gas_so] -= g_dso;
00579                     pp_block_p[p3_gas_sg] -= g_dsg;
00580                     ps_block_p[p3_gas] += -g_dsw;
00581 
00582                     pp_block_m_up[p3_gas_po] += g_dpo_up;
00583                     pp_block_m_up[p3_gas_so] += g_dso_up;
00584                     pp_block_m_up[p3_gas_sg] += g_dsg_up;
00585                     ps_block_m_up[p3_gas] += g_dsw_up;
00586 
00587                     pp_block_p_up[p3_gas_po] -= g_dpo_up;
00588                     pp_block_p_up[p3_gas_so] -= g_dso_up;
00589                     pp_block_p_up[p3_gas_sg] -= g_dsg_up;
00590                     ps_block_p_up[p3_gas] += -g_dsw_up;
00591                   }
00592                 else if (is_w)
00593                   {
00594                     pp_block_m[p2ow_oil_po] += dpo;
00595                     pp_block_m[p2ow_oil_so] += dso;
00596                     ps_block_m[p2ow_oil] += dsw;
00597 
00598                     pp_block_p[p2ow_oil_po] -= dpo;
00599                     pp_block_p[p2ow_oil_so] -= dso;
00600                     ps_block_p[p2ow_oil] += -dsw;
00601 
00602                     pp_block_m_up[p2ow_oil_po] += dpo_up;
00603                     pp_block_m_up[p2ow_oil_so] += dso_up;
00604                     ps_block_m_up[p2ow_oil] += dsw_up;
00605 
00606                     pp_block_p_up[p2ow_oil_po] -= dpo_up;
00607                     pp_block_p_up[p2ow_oil_so] -= dso_up;
00608                     ps_block_p_up[p2ow_oil] += -dsw_up;
00609                   }
00610                 else if (is_g)
00611                   {
00612                     // oil
00613                     pp_block_m[p2og_oil_po] += dpo;
00614                     pp_block_m[p2og_oil_sg] += dsg;
00615                     ps_block_m[p2og_oil] += dso;
00616 
00617                     pp_block_p[p2og_oil_po] -= dpo;
00618                     pp_block_p[p2og_oil_sg] -= dsg;
00619                     ps_block_p[p2og_oil] += -dso;
00620 
00621                     pp_block_m_up[p2og_oil_po] += dpo_up;
00622                     pp_block_m_up[p2og_oil_sg] += dsg_up;
00623                     ps_block_m_up[p2og_oil] += dso_up;
00624 
00625                     pp_block_p_up[p2og_oil_po] -= dpo_up;
00626                     pp_block_p_up[p2og_oil_sg] -= dsg_up;
00627                     ps_block_p_up[p2og_oil] += -dso_up;
00628 
00629                     // gas
00630                     pp_block_m[p2og_gas_po] += g_dpo;
00631                     pp_block_m[p2og_gas_sg] += g_dsg;
00632                     ps_block_m[p2og_gas] += g_dso;
00633 
00634                     pp_block_p[p2og_gas_po] -= g_dpo;
00635                     pp_block_p[p2og_gas_sg] -= g_dsg;
00636                     ps_block_p[p2og_gas] += -g_dso;
00637 
00638                     pp_block_m_up[p2og_gas_po] += g_dpo_up;
00639                     pp_block_m_up[p2og_gas_sg] += g_dsg_up;
00640                     ps_block_m_up[p2og_gas] += g_dso_up;
00641 
00642                     pp_block_p_up[p2og_gas_po] -= g_dpo_up;
00643                     pp_block_p_up[p2og_gas_sg] -= g_dsg_up;
00644                     ps_block_p_up[p2og_gas] += -g_dso_up;
00645                   }
00646                 else
00647                   {
00648                     pp_block_m[0] += dpo;
00649                     pp_block_p[0] -= dpo;
00650 
00651                     pp_block_m_up[0] += dpo_up;
00652                     pp_block_p_up[0] -= dpo_up;
00653                   }
00654                 if (n_phases > 1)
00655                   {
00656                     M_MINUS_VV_PROD (n_phases, ps_block_m_up, sp_block_up, pp_block_m_up);
00657                     V_MINUS_VS_PROD (n_phases, ps_block_m_up, s_rhs_block_up, rhs_m);
00658 
00659                     M_MINUS_VV_PROD (n_phases, ps_block_p_up, sp_block_up, pp_block_p_up);
00660                     V_MINUS_VS_PROD (n_phases, ps_block_p_up, s_rhs_block_up, rhs_p);
00661                   }
00662               }
00663             // gas
00664             if (is_g)
00665               {
00666                 //index_t mobg = 0;
00667                 k = k_j + d_g;
00668 
00669                 ps_block_p_up[0] = ps_block_p_up[1] = ps_block_p_up[2] = 0;
00670                 ps_block_m_up[0] = ps_block_m_up[1] = ps_block_m_up[2] = 0;
00671                 sp_block_up = sp_diag_val + sp_size_g;
00672                 s_rhs_block_up = s_rhs_val + rhs_size_g;
00673 
00674                 if (up_cell_dg == cell_ind_block[0])
00675                   {
00676                     pp_block_m_up = &m_mem[0];
00677                     pp_block_p_up = &p_mem[0];
00678                   }
00679                 else
00680                   {
00681                     pp_block_m_up = &m_mem[9];
00682                     pp_block_p_up = &p_mem[9];
00683                   }
00684 
00685                 // calculate pressure derivaties
00686                 dpo = truns_mob_g * truns_block[j];
00687                 dpo_up = 0;
00688                 if (ic == up_cell_dg)
00689                   dpo += truns_p_mob_g * truns_block[j] * cell_pot[k];
00690                 else
00691                   dpo_up = truns_p_mob_g * truns_block[j] * cell_pot[k];
00692 
00693                 dsw = 0;
00694                 dso = 0;
00695                 dsg = 0;
00696                 dsw_up = dso_up = dsg_up = 0;
00697                 // calculate mobility indexies
00698                 //mobg  = d_g * n_phases + d_g;
00699 
00700                 if (is_o)
00701                   {
00702                     dsg += truns_mob_g * truns_block[j] * data_[ic].s_deriv_cap_pressure[d_g];
00703                     dsg_up = 0;
00704                     if (ic == up_cell_dg)
00705                       {
00706                         dsg += truns_s_mob_gg * truns_block[j] * cell_pot[k];
00707                       }
00708                     else
00709                       {
00710                         dsg_up = truns_s_mob_gg * truns_block[j] * cell_pot[k];
00711                       }
00712                   }
00713                 if (n_phases > 1 && !(FI_CHK_SG (main_vars_, up_cell_dg)))
00714                   {
00715                     dpo = dpo_up = 0;
00716                     dso = dso_up = 0;
00717                     dsg = dsg_up = 0;
00718                     dsw = dsw_up = 0;
00719                   }
00720 
00721                 if (n_phases == 3)
00722                   {
00723                     pp_block_m[p3_gas_po] += dpo;
00724                     pp_block_m[p3_gas_so] += dso;
00725                     pp_block_m[p3_gas_sg] += dsg;
00726                     ps_block_m[p3_gas] += dsw;
00727 
00728                     pp_block_p[p3_gas_po] -= dpo;
00729                     pp_block_p[p3_gas_so] -= dso;
00730                     pp_block_p[p3_gas_sg] -= dsg;
00731                     ps_block_p[p3_gas] += -dsw;
00732 
00733                     // up block
00734                     pp_block_m_up[p3_gas_po] += dpo_up;
00735                     pp_block_m_up[p3_gas_so] += dso_up;
00736                     pp_block_m_up[p3_gas_sg] += dsg_up;
00737                     ps_block_m_up[p3_gas] += dsw_up;
00738 
00739                     pp_block_p_up[p3_gas_po] -= dpo_up;
00740                     pp_block_p_up[p3_gas_so] -= dso_up;
00741                     pp_block_p_up[p3_gas_sg] -= dsg_up;
00742                     ps_block_p_up[p3_gas] += -dsw_up;
00743                   }
00744                 else if (n_phases == 2)
00745                   {
00746                     pp_block_m[p2og_gas_po] += dpo;
00747                     pp_block_m[p2og_gas_sg] += dsg;
00748                     ps_block_m[p2og_gas] += dso;
00749 
00750                     pp_block_p[p2og_gas_po] -= dpo;
00751                     pp_block_p[p2og_gas_sg] -= dsg;
00752                     ps_block_p[p2og_gas] += -dso;
00753 
00754                     // up block
00755                     pp_block_m_up[p2og_gas_po] += dpo_up;
00756                     pp_block_m_up[p2og_gas_sg] += dsg_up;
00757                     ps_block_m_up[p2og_gas] += dso_up;
00758 
00759                     pp_block_p_up[p2og_gas_po] -= dpo_up;
00760                     pp_block_p_up[p2og_gas_sg] -= dsg_up;
00761                     ps_block_p_up[p2og_gas] += -dso_up;
00762                   }
00763                 else
00764                   {
00765                     pp_block_m[0] += dpo;
00766                     pp_block_p[0] -= dpo;
00767 
00768                     // up block
00769                     pp_block_m_up[0] += dpo_up;
00770                     pp_block_p_up[0] -= dpo_up;
00771                   }
00772                 if (n_phases > 1)
00773                   {
00774                     M_MINUS_VV_PROD (n_phases, ps_block_m_up, sp_block_up, pp_block_m_up);
00775                     M_MINUS_VV_PROD (n_phases, ps_block_p_up, sp_block_up, pp_block_p_up);
00776                     V_MINUS_VS_PROD (n_phases, ps_block_m_up, s_rhs_block_up, rhs_m);
00777                     V_MINUS_VS_PROD (n_phases, ps_block_p_up, s_rhs_block_up, rhs_p);
00778                   }
00779               }
00780             M_MINUS_VV_PROD (n_phases, ps_block_m, sp_block, pp_block_m);
00781             M_MINUS_VV_PROD (n_phases, ps_block_p, sp_block, pp_block_p);
00782 
00783             V_MINUS_VS_PROD (n_phases, ps_block_m, s_rhs_block, rhs_m);
00784             V_MINUS_VS_PROD (n_phases, ps_block_p, s_rhs_block, rhs_p);
00785           }
00786 
00787       }
00788 
00789     public:
00790 
00791       const fi_operator_impl <strategy_t> &fi_operator_;
00792       index_t                               n_phases;
00793       bool                                  is_w;
00794       bool                                  is_g;
00795       bool                                  is_o;
00796       index_t                               d_w;
00797       index_t                               d_g;
00798       index_t                               d_o;
00799       const item_array_t                    &saturation_3p_;
00800       const item_array_t                    &pressure_;
00801       const item_array_t                    &gas_oil_ratio_;
00802       const main_var_array_t                &main_vars_;
00803       const data_array_t                    &data_;
00804       item_t                                gravity_;
00805       index_t                               n_sec_vars;
00806       item_array_t                          &rhs_;
00807       item_array_t                          &reg_values_;
00808     };
00809   } // tpfa
00810 
00811   template <typename strategy_t>
00812   bool
00813   fi_operator_impl <strategy_t>::block_connections_mpfa (const item_t &dt)
00814   {
00817     //if (n_phases == 3)
00818     //  {
00819     //    return mpfa_calc_3phase <strategy_t> (this).calc (dt);
00820     //  }
00821     //else if (n_phases == 2 && is_w)
00822     //  {
00823     //    return mpfa_calc_2phase_ow <strategy_t> (this).calc (dt);
00824     //  }
00825     //else if (n_phases == 2 && is_g)
00826     //  {
00827     //    return mpfa_calc_2phase_og <strategy_t> (this).calc (dt);
00828     //  }
00829     //else if (is_w)
00830     //  {
00831     //    return mpfa_calc_1phase_w <strategy_t> (this).calc (dt);
00832     //  }
00833     //else if (is_g)
00834     //  {
00835     //    return mpfa_calc_1phase_g <strategy_t> (this).calc (dt);
00836     //  }
00837     //else if (is_o)
00838     //  {
00839     //    return mpfa_calc_1phase_o <strategy_t> (this).calc (dt);
00840     //  }
00841     //else
00842     //  {
00843     //    throw bs_exception ("block_connections_mpfa", "NOT IMPL YET");
00844     //  }
00845 
00846     index_t b_sqr = n_phases * n_phases;
00847     index_t equ_w = 0, equ_g = 0, equ_o = 0;
00848 
00849     boost::array <item_t, FI_PHASE_TOT> ave_density;     // average density of phase over all cells involved in MPFA connection
00850 
00851     //potential
00852     item_t cell_pot[2 * FI_PHASE_TOT];    // potential for cell
00853     item_t sum_cell_pot[FI_PHASE_TOT];
00854 
00855     //up_stream
00856     typedef tpfa::mpfa_impl <strategy_t> mpfa_impl_t;
00857     typename mpfa_impl_t::up_cell_array_t up_cell;    // index of upstream cell
00858 
00859     //rhs
00860     item_t * rhs_m;
00861     item_t * rhs_p;
00862 
00863     //jacobian
00864     boost::array <item_t, 2 * 9> m_mem; 
00865     boost::array <item_t, 2 * 9> p_mem; 
00866 
00867     assign (flux_rhs_, 0);
00868 
00869     // TPFA: (SHOULD BE REMOVED AFTER FULL MPFA IMPL)
00870     // number of cells used to build MPFA for one connection
00871     index_t n_cells_in_conn = 2;          
00872 
00873     if (n_phases == 3)
00874       {
00875         equ_w = p3_wat;
00876         equ_g = p3_gas;
00877         equ_o = p3_oil;
00878       }
00879     else if (n_phases == 2 && is_w)
00880       {
00881         equ_w = p2ow_wat;
00882         equ_o = p2ow_oil;
00883       }
00884     else if (n_phases == 2 && is_g)
00885       {
00886         equ_g = p2og_gas;
00887         equ_o = p2og_oil;
00888       }
00889     else
00890       {
00891         equ_g = equ_w = equ_o = 0;
00892       }
00893 
00894     mpfa_impl_t mpfa (*this, flux_rhs_, reg_values_);
00895 
00896     // main loop through connections
00897     for (index_t i = 0; i < n_connections_; ++i)
00898       {
00899         index_t row_i = trns_rows_ptr_[i];
00900         const item_t *truns_block = &trns_values_[row_i];       // temporary array for TPFA connection (SHOULD BE REMOVED AFTER FULL MPFA IMPL)
00901 
00902         // temporary array for TPFA cell ind (SHOULD BE REMOVED AFTER FULL MPFA IMPL)
00903         const index_t *cell_ind_block = &trns_cols_ptr_[row_i];
00904 
00905         index_t cell_m = cell_ind_block[0];                           //pervii v cols index v stroke          // index of primary cell in MPFA connection (minus)
00906         index_t cell_p = cell_ind_block[1];                           //vtoroi          // index of primary cell in MPFA connection (plus)
00907         const index_t *m_mem_ptr = &m_array_[row_i];      // ptr to the block in matrix row cell_m for columns with indexies #cell_ind_block
00908         const index_t *p_mem_ptr = &p_array_[row_i];      // ptr to the block in matrix row cell_p for columns with indexies #cell_ind_block
00909 
00910         //fills ave_density and ++s counter
00911         mpfa.mpfa_calc_avarage_density(&ave_density[0], cell_ind_block, n_cells_in_conn);
00912 
00913         //fills cell_pot and sum_cell_pot
00914         mpfa.mpfa_calc_potential(cell_pot, sum_cell_pot, &ave_density[0], truns_block, depths_,
00915                             cell_ind_block, cell_m, dt, n_cells_in_conn);
00916 
00917         //fills up_cell
00918         mpfa.mpfa_calc_upstream(up_cell, sum_cell_pot, cell_m, cell_p);
00919 
00920         //fills rhs_m and rhs_p
00921         mpfa.mpfa_fill_rhs(rhs_m, rhs_p, up_cell, truns_block, cell_pot, mpfa.rhs_, dt, n_cells_in_conn,
00922                       cell_m, cell_p, equ_w, equ_g, equ_o);
00923 
00924 
00925         //fills m_mem and p_mem
00926         mpfa.mpfa_fill_jacobian(m_mem, p_mem, up_cell, cell_pot, truns_block, rhs_m, rhs_p,
00927                            cell_ind_block, sp_diag_, s_rhs_, dt, n_cells_in_conn,
00928                            cell_m, cell_p, equ_w, equ_g, equ_o);
00929 
00930         item_t * block_m1 = &mpfa.reg_values_[m_mem_ptr[0] * b_sqr];
00931         item_t * block_p1 = &mpfa.reg_values_[p_mem_ptr[0] * b_sqr];
00932         item_t * block_m2 = &mpfa.reg_values_[m_mem_ptr[1] * b_sqr];
00933         item_t * block_p2 = &mpfa.reg_values_[p_mem_ptr[1] * b_sqr];
00934 
00935         for (index_t j = 0; j < b_sqr; ++j)
00936           {
00937             block_m1[j] += m_mem[j];
00938             block_p1[j] += p_mem[j];
00939 
00940             block_m2[j] += m_mem[9 + j];
00941             block_p2[j] += p_mem[9 + j];
00942           }
00943 
00944         // END OF JACOBIAN filling
00945       }
00946 
00947     return 0;
00948   }
00949 
00950 } // namespace blue_sky
00951 
00952 
00953 #endif  // BS_FI_OPERATOR_BLOCK_CONNECTIONS_MPFA_H_
00954
fi_operator_block_connections_mpfa.h
