From 6fc24236dff0d41620deb23c6fc2d350ffd44e4f Mon Sep 17 00:00:00 2001 From: Liu Ming Date: Mon, 9 Dec 2013 22:57:44 +0800 Subject: [PATCH] fix FullyImplicitTwophasePolymersolver constructor problem. add some input commits for debugging. --- examples/.sim_poly2p_fincomp_ad.cpp.swp | Bin 0 -> 20480 bytes examples/sim_2p_fim.cpp | 262 ++++++++++++++++++ examples/sim_2p_fincomp_ad.cpp | 91 ++++++ examples/sim_poly2p_fincomp_ad.cpp | 146 ++++++++++ .../FullyImplicitTwophasePolymerSolver.cpp | 80 ++++-- .../FullyImplicitTwophasePolymerSolver.hpp | 8 +- opm/polymer/fullyimplicit/PolymerPropsAd.hpp | 1 - 7 files changed, 561 insertions(+), 27 deletions(-) create mode 100644 examples/.sim_poly2p_fincomp_ad.cpp.swp create mode 100644 examples/sim_2p_fim.cpp create mode 100644 examples/sim_2p_fincomp_ad.cpp create mode 100644 examples/sim_poly2p_fincomp_ad.cpp diff --git a/examples/.sim_poly2p_fincomp_ad.cpp.swp b/examples/.sim_poly2p_fincomp_ad.cpp.swp new file mode 100644 index 0000000000000000000000000000000000000000..00e8add3fc22b334f8bb9566a5eba808f776ff56 GIT binary patch literal 20480 zcmeHPYm6jS6)wORFC`KJ{tz#{WT~0yt?r(6VdD6u0YpQ5{XnDNxsR@@?&_IE zK#8eNzNzke&pqedb02ljxo3B|ee1$ecAfQJ1D{I_q9T5&%l?CHnBR>+eU-*m!$z|&Uhhc*??_U)b@xpCOH z9nW(6{T&6=1hoZf3shTRkoETLzu4G2+i4TQuEtgD(s!P!rmahB3)B{fE<^=X><|LUGu(w)(kdP+sQZxu{mj^+*r*k_y!1i>kxnW(KN`?iQ*X^!z zH;a>bFCHuh{2&7Zk@A5&5Oc_AN-?OW&*_?sx;~jW4*7~BtCcAhlenKA@`h&U({Ci1 znA$F)K(Qnfj(l8u-3d~A-A~;nWAv~?jEX32GF9t%X_MiPD4om;M$UZKcAYev7o%6u zktdnY({A^Wm&Qpy^W*3SzUr*|@gO;95^F#5_~w)x=92DK$}v4cep4vT!#)bNTrm#p_nnGkZsZ-f>(zv>>OS3fG7?bww zzL6R+HK?XDS|g?iY9Zi#Q?dMPN2$SowFuoVK##&9}=cgHHX z%%^>2)4W=?WpgyBc_|uo*owc-Qx?U{U3H=r&isgBp@iI7_5&<6gQkK~w6y3dC6*{4 zv~6|vji_8$#+5Acbxx9LRAO54OBl$pf&*+Hn`SPLG9S;y*iX4!gr#P*peN9!6&v)h zt=Ojft4_+X8+NSrv=B5c#HdbHwO4yoxpiB-f~0bT6q$)EazdW=9hb`)(@%KM-`s|P z-i)%ECMCUTq-rqYY`I*#%22$L_?`vAWuDk^Prl(Kr5i_S#!lf+GH|m&!o4G)-7Q+6 z;WV|xE*d#2Jh5DxewwPsM0ixz2y#tkM%b<+?99=f*q!m+qh0#W^5A=Sf;$F7YN2P< z7Ii_dU0t|YEEGMaww*#xsoCgyG#k+trroUGb3{X8RYAB}Bs}Z}em~{ZyF~{%p$XnF z?iw17RhU$P?9KprqaiQO16Y#=wp3pfJuha9q{60AqH8EdcN=1sm_kmR^_9S7UK7qdDS50;@pCFs=GGD zP7A@yaU8Hx&9DPl2NQG_w&fX5x?R`~M#0x3Nj&IFtD8=;g0lwxt~Xf^hYxMEighq& zoHnv#%b;-!9c)D=eUl8ZKPA)Cg3pw(3)!{7z~gKldAKbX2Z5F!x40a|*4Us~tNs4L zsz!rM>~gJPAJSI24?DIzOapS%ag#XRk@ghU@Xk~r_PkRXE8oeft>c_MlAM?0{sYcbr@@n9p_+yep^CowYaw64q*Sxibc|8Dwg!H4RR+e7WpyWIDToF`tcZ91MjY$zhbR6 z4D$aUfo~kcZzlhLM#Uf=hJXKQU=_Fo_!E5l$AB%M3tR&H8h-t^fiu8Ez`ei#_z-Y4 z@G|`SzXIO?9t0i$Qs5-e1}+1hgOC41;4a{Uz#iaf`1qd%2Ebw9UBC;_?+?K5fbRfz z0p#OZSnJsUwkFl+QxGB)Y3QH;Ri@IdNi*Izy{w78z+8AF@ZSsf90i?~ zyQXeq=$c2mfVVoqfTyNAND_+9h2WBZpAo$Cdy@jcoQTaFV4jM2_2NWE0(=Z&dyb)J zA_68NBBX6P$7W`H9mUXb)uC|F$`C!IfFF4kh}1~mqG%C4uW6ECR#LoA&Iioe1I*Rc z?hmkZJb0X_aPj*AW3i#8UH_L_;)oWr;{B zCK*N1b!1UC*U`$XvXw_SAf2`~h2S)-_%T!{nbXe-1bT{0n~leZVcii}2^q0QUpr;~xgD2Ix-ztH5dC zW56NcYT(cC<(~z<349SqflGm}BSygi8#oCt;2hwmxWoSva3637-~;CaqCaoYLqM%o zZGqYXwFPPmydn#1{{qyt6JCU+RH@fHQ>u!;K?+tUO0iacvLOnY)zP*}P zWL9;^0ew#QI%iGog2DZs=yZNeB@>K%RCW#J3%u%7D)Jj&qhu=5vV3HEZ4;`9|Dw#x zgU`QBQWX*E%%}0hs#5J(a#g9eoWny^RZ~QZ8cTJe26-mbNTSt9WjWbur1I?a$VF9Z zLenKrx%+P)?UTBCHTHAt|fSx|CxzzOCog(a?XZi9OdsX{B5&>l&oY_cBl4J2mcbFN~bU@ea90f zgy2-XM0HME!ULQ5+Py)9SJ*JQmR6q#K*$W*6T`45?f0ojW9W8|^I#EeZ$eXK92FR0 s(5hR=7=nxns8v?%71|+I(w+_x!}O;;CvjE8MANFM2 + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include + +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include + + +namespace +{ + void warnIfUnusedParams(const Opm::parameter::ParameterGroup& param) + { + if (param.anyUnused()) { + std::cout << "-------------------- Unused parameters: --------------------\n"; + param.displayUsage(); + std::cout << "----------------------------------------------------------------" << std::endl; + } + } +} // anon namespace + + + +// ----------------- Main program ----------------- +int +main(int argc, char** argv) +try +{ + using namespace Opm; + + std::cout << "\n================ Test program for incompressible two-phase flow ===============\n\n"; + parameter::ParameterGroup param(argc, argv, false); + std::cout << "--------------- Reading parameters ---------------" << std::endl; + + + // If we have a "deck_filename", grid and props will be read from that. + bool use_deck = param.has("deck_filename"); + boost::scoped_ptr deck; + boost::scoped_ptr grid; + boost::scoped_ptr props; + TwophaseState state; + double gravity[3] = { 0.0 }; + if (use_deck) { + std::string deck_filename = param.get("deck_filename"); + deck.reset(new EclipseGridParser(deck_filename)); + // Grid init + grid.reset(new GridManager(*deck)); + // Rock and fluid init + props.reset(new IncompPropsAdFromDeck(*deck, *grid->c_grid())); + // Gravity. + gravity[2] = deck->hasField("NOGRAV") ? 0.0 : unit::gravity; + // Init state variables (saturation and pressure). + int num_cells = grid->c_grid()->number_of_cells; + if (param.has("init_saturation")) { + //initStateBasic(*grid->c_grid(), *props, param, gravity[2], state); + const double init_saturation = param.get("init_saturation"); + for (int c = 0; c < num_cells; ++c) { + state.saturation()[2*c] = init_saturation; + state.saturation()[2*c+1] = 1. - init_saturation; + } + } else { + if (deck->hasField("PRESSURE")) { + // Set saturations from SWAT/SGAS, pressure from PRESSURE. + std::vector& s = state.saturation(); + std::vector& p = state.pressure(); + const std::vector& p_deck = deck->getFloatingPointValue("PRESSURE"); + // water-oil or water-gas: we require SWAT + if (!deck->hasField("SWAT")) { + OPM_THROW(std::runtime_error, "initStateFromDeck(): missing SWAT keyword in 2-phase init"); + } + const std::vector& sw_deck = deck->getFloatingPointValue("SWAT"); + for (int c = 0; c < num_cells; ++c) { + int c_deck = (grid->c_grid()->global_cell == NULL) ? c : grid->c_grid()->global_cell[c]; + s[2*c] = sw_deck[c_deck]; + s[2*c + 1] = 1.0 - sw_deck[c_deck]; + p[c] = p_deck[c_deck]; + } + } + } + } else { + // Grid init. + const int nx = param.getDefault("nx", 100); + const int ny = param.getDefault("ny", 100); + const int nz = param.getDefault("nz", 1); + const double dx = param.getDefault("dx", 1.0); + const double dy = param.getDefault("dy", 1.0); + const double dz = param.getDefault("dz", 1.0); + grid.reset(new GridManager(nx, ny, nz, dx, dy, dz)); + // Rock and fluid init. + props.reset(new IncompPropsAdBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells)); + // Rock compressibility. + // Gravity. + gravity[2] = param.getDefault("gravity", 0.0); + int num_cells = grid->c_grid()->number_of_cells; + } + + // Warn if gravity but no density difference. + bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0); + const double *grav = use_gravity ? &gravity[0] : 0; + + // Initialising src + std::vector src(num_cells, 0.0); + if (use_deck) { + // Do nothing, wells will be the driving force, not source terms. + if (deck->hasField("SRC")) { + const std::vector& src_deck = deck->getFloatingPointValue("SRC"); + for (int c = 0; c < num_cells; ++c) { + int c_deck = (grid->c_grid()->global_cell == NULL) ? c : grid->c_grid()->global_cell[c]; + src[c] = src_deck[c_deck]; + } + } + } else { + // Compute pore volumes, in order to enable specifying injection rate + // terms of total pore volume. + std::vector porevol; + computePorevolume(*grid->c_grid(), props->porosity(), porevol); + const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0); + const double default_injection = use_gravity ? 0.0 : 0.1; + const double flow_per_sec = param.getDefault("injected_porevolumes_per_day", default_injection) + *tot_porevol_init/unit::day; + src[0] = flow_per_sec; + src[num_cells - 1] = -flow_per_sec; + } + + in: num_cells = grid->c_grid()->number_of_cells; + + // Linear solver. + LinearSolverFactory linsolver(param); + + // Write parameters used for later reference. + bool output = param.getDefault("output", true); + std::ofstream epoch_os; + std::string output_dir; + if (output) { + output_dir = + param.getDefault("output_dir", std::string("output")); + boost::filesystem::path fpath(output_dir); + try { + create_directories(fpath); + } + catch (...) { + OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath); + } + std::string filename = output_dir + "/epoch_timing.param"; + epoch_os.open(filename.c_str(), std::fstream::trunc | std::fstream::out); + // open file to clean it. The file is appended to in SimulatorTwophase + filename = output_dir + "/step_timing.param"; + std::fstream step_os(filename.c_str(), std::fstream::trunc | std::fstream::out); + step_os.close(); + param.writeParam(output_dir + "/simulation.param"); + } + + + std::cout << "\n\n================ Starting main simulation loop ===============\n" + << " (number of epochs: " + << (use_deck ? deck->numberOfEpochs() : 1) << ")\n\n" << std::flush; + + SimulatorReport rep; + if (!use_deck) { + // Simple simulation without a deck. + SimulatorFullyImplicitTwophase simulator(param, + *grid->c_grid(), + *props, + linsolver, + src); + SimulatorTimer simtimer; + simtimer.init(param); + warnIfUnusedParams(param); + rep = simulator.run(simtimer, state, src); + } else { + // With a deck, we may have more epochs etc. + int step = 0; + SimulatorTimer simtimer; + // Use timer for last epoch to obtain total time. + deck->setCurrentEpoch(deck->numberOfEpochs() - 1); + simtimer.init(*deck); + const double total_time = simtimer.totalTime(); + for (int epoch = 0; epoch < deck->numberOfEpochs(); ++epoch) { + // Set epoch index. + deck->setCurrentEpoch(epoch); + + // Update the timer. + if (deck->hasField("TSTEP")) { + simtimer.init(*deck); + } else { + if (epoch != 0) { + OPM_THROW(std::runtime_error, "No TSTEP in deck for epoch " << epoch); + } + simtimer.init(param); + } + simtimer.setCurrentStepNum(step); + simtimer.setTotalTime(total_time); + + // Report on start of epoch. + std::cout << "\n\n-------------- Starting epoch " << epoch << " --------------" + << "\n (number of steps: " + << simtimer.numSteps() - step << ")\n\n" << std::flush; + + + // Create and run simulator. + SimulatorFullyImplicitTwophase simulator(param, + *grid->c_grid(), + *props, + linsolver, + src); + if (epoch == 0) { + warnIfUnusedParams(param); + } + SimulatorReport epoch_rep = simulator.run(simtimer, state, src); + if (output) { + epoch_rep.reportParam(epoch_os); + } + // Update total timing report and remember step number. + rep += epoch_rep; + step = simtimer.currentStepNum(); + } + } + + std::cout << "\n\n================ End of simulation ===============\n\n"; + rep.report(std::cout); + + if (output) { + std::string filename = output_dir + "/walltime.param"; + std::fstream tot_os(filename.c_str(),std::fstream::trunc | std::fstream::out); + rep.reportParam(tot_os); + } + +} +catch (const std::exception &e) { + std::cerr << "Program threw an exception: " << e.what() << "\n"; + throw; +} + diff --git a/examples/sim_2p_fincomp_ad.cpp b/examples/sim_2p_fincomp_ad.cpp new file mode 100644 index 000000000..544ccb8c2 --- /dev/null +++ b/examples/sim_2p_fincomp_ad.cpp @@ -0,0 +1,91 @@ +#include "config.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +#include +#include + +#include +#include +#include + +int main (int argc, char** argv) +try +{ + int nx = 20; + int ny = 20; + int nz = 1; + double dx = 10.0; + double dy = 10.0; + double dz = 10.0; + using namespace Opm; + parameter::ParameterGroup param(argc, argv, false); + GridManager grid_manager(nx, ny, nz, dx, dy, dz); + const UnstructuredGrid& grid = *grid_manager.c_grid(); + int num_cells = grid.number_of_cells; + int num_phases = 2; + using namespace Opm::unit; + using namespace Opm::prefix; + std::vector density(num_phases, 1000.0); + std::vector viscosity(num_phases, 1.0*centi*Poise); + double porosity = 0.5; + double permeability = 10.0*milli*darcy; + SaturationPropsBasic::RelPermFunc rel_perm_func = SaturationPropsBasic::Linear; + IncompPropsAdBasic props(num_phases, rel_perm_func, density, viscosity, + porosity, permeability, grid.dimensions, num_cells); + std::vector omega; + std::vector src(num_cells, 0.0); + src[0] = 1.; + src[num_cells-1] = -1.; + + FlowBCManager bcs; + LinearSolverUmfpack linsolver; + FullyImplicitTwoPhaseSolver solver(grid, props, linsolver); + std::vector porevol; + Opm::computePorevolume(grid, props.porosity(), porevol); + const double dt = param.getDefault("dt", 0.1) * day; + const int num_time_steps = param.getDefault("nsteps", 20); + std::vector allcells(num_cells); + for (int cell = 0; cell < num_cells; ++cell) { + allcells[cell] = cell; + } + TwophaseState state; + state.init(grid, 2); + + //initial sat + for (int c = 0; c < num_cells; ++c) { + state.saturation()[2*c] = 0.2; + state.saturation()[2*c+1] = 0.8; + } + std::vector p(num_cells, 200*Opm::unit::barsa); + state.pressure() = p; +// state.setFirstSat(allcells, props, TwophaseState::MinSat); + std::ostringstream vtkfilename; + + for (int i = 0; i < num_time_steps; ++i) { + solver.step(dt, state, src); + vtkfilename.str(""); + vtkfilename << "sim_2p_fincomp-" << std::setw(3) << std::setfill('0') << i << ".vtu"; + std::ofstream vtkfile(vtkfilename.str().c_str()); + Opm::DataMap dm; + dm["saturation"] = &state.saturation(); + dm["pressure"] = &state.pressure(); + Opm::writeVtkData(grid, dm, vtkfile); + } +} +catch (const std::exception &e) { + std::cerr << "Program threw an exception: " << e.what() << "\n"; + throw; +} diff --git a/examples/sim_poly2p_fincomp_ad.cpp b/examples/sim_poly2p_fincomp_ad.cpp new file mode 100644 index 000000000..80e5f3e6a --- /dev/null +++ b/examples/sim_poly2p_fincomp_ad.cpp @@ -0,0 +1,146 @@ +#include "config.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +int main (int argc, char** argv) +try +{ + using namespace Opm; + parameter::ParameterGroup param(argc, argv, false); + bool use_poly_deck = param.has("deck_filename"); + if (!use_poly_deck) { + OPM_THROW(std::runtime_error, "Polymer Properties must be read from deck_filename\n"); + } + std::string deck_filename = param.get("deck_filename"); + EclipseGridParser deck = EclipseGridParser(deck_filename); + int nx = param.getDefault("nx", 20); + int ny = param.getDefault("ny", 20); + int nz = 1; + double dx = 2.0; + double dy = 2.0; + double dz = 0.5; + GridManager grid_manager(nx, ny, nz, dx, dy, dz); + const UnstructuredGrid& grid = *grid_manager.c_grid(); + int num_cells = grid.number_of_cells; + int num_phases = 2; + using namespace Opm::unit; + using namespace Opm::prefix; + std::vector density(num_phases, 1000.0); + std::vector viscosity(num_phases, 1.0*centi*Poise); + viscosity[0] = 0.5 * centi * Poise; + viscosity[0] = 5 * centi * Poise; + double porosity = 0.35; + double permeability = 10.0*milli*darcy; + SaturationPropsBasic::RelPermFunc rel_perm_func = SaturationPropsBasic::Linear; + IncompPropsAdBasic props(num_phases, rel_perm_func, density, viscosity, + porosity, permeability, grid.dimensions, num_cells); + + // Init polymer properties. + // Setting defaults to provide a simple example case. + PolymerProperties polymer_props(deck); + #if 0 + if (use_poly_deck) { + } else { + double c_max = param.getDefault("c_max_limit", 5.0); + double mix_param = param.getDefault("mix_param", 1.0); + double rock_density = param.getDefault("rock_density", 1000.0); + double dead_pore_vol = param.getDefault("dead_pore_vol", 0.15); + double res_factor = param.getDefault("res_factor", 1.) ; // res_factor = 1 gives no change in permeability + double c_max_ads = param.getDefault("c_max_ads", 1.); + int ads_index = param.getDefault("ads_index", Opm::PolymerProperties::NoDesorption); + std::vector c_vals_visc(2, -1e100); + c_vals_visc[0] = 0.0; + c_vals_visc[1] = 7.0; + std::vector visc_mult_vals(2, -1e100); + visc_mult_vals[0] = 1.0; + // poly_props.visc_mult_vals[1] = param.getDefault("c_max_viscmult", 30.0); + visc_mult_vals[1] = 20.0; + std::vector c_vals_ads(3, -1e100); + c_vals_ads[0] = 0.0; + c_vals_ads[1] = 2.0; + c_vals_ads[2] = 8.0; + std::vector ads_vals(3, -1e100); + ads_vals[0] = 0.0; + ads_vals[1] = 0.0015; + ads_vals[2] = 0.0025; + PolymerProperties polymer_props; + polymer_props.set(c_max, mix_param, rock_density, dead_pore_vol, res_factor, c_max_ads, + static_cast(ads_index), + c_vals_visc, visc_mult_vals, c_vals_ads, ads_vals); + } + #endif + PolymerPropsAd polymer_props_ad(polymer_props); + std::vector omega; + std::vector src(num_cells, 0.0); + std::vector src_polymer(num_cells); + src[0] = 10. / day; + src[num_cells-1] = -1. / day; + + PolymerInflowBasic polymer_inflow(param.getDefault("poly_start_days", 30.0)*Opm::unit::day, + param.getDefault("poly_end_days", 80.0)*Opm::unit::day, + param.getDefault("poly_amount", polymer_props.cMax())); + FlowBCManager bcs; + LinearSolverUmfpack linsolver; + FullyImplicitTwophasePolymerSolver solver(grid, props,polymer_props_ad, linsolver); + std::vector porevol; + Opm::computePorevolume(grid, props.porosity(), porevol); + const double dt = param.getDefault("dt", 10.) * day; + const int num_time_steps = param.getDefault("nsteps", 10); + std::vector allcells(num_cells); + for (int cell = 0; cell < num_cells; ++cell) { + allcells[cell] = cell; + } + PolymerState state; + state.init(grid, 2); + //initial sat + for (int c = 0; c < num_cells; ++c) { + state.saturation()[2*c] = 0.2; + state.saturation()[2*c+1] = 0.8; + } + std::vector p(num_cells, 200*Opm::unit::barsa); + state.pressure() = p; + + std::vector c(num_cells, 0.0); + state.concentration() = c; + std::ostringstream vtkfilename; + double currentime = 0; + for (int i = 0; i < num_time_steps; ++i) { + currentime += dt; + polymer_inflow.getInflowValues(currentime, currentime+dt, src_polymer); + solver.step(dt, state, src, src_polymer); + vtkfilename.str(""); + vtkfilename << "sim_poly2p_fincomp_ad_" << std::setw(3) << std::setfill('0') << i << ".vtu"; + std::ofstream vtkfile(vtkfilename.str().c_str()); + Opm::DataMap dm; + dm["saturation"] = &state.saturation(); + dm["pressure"] = &state.pressure(); + Opm::writeVtkData(grid, dm, vtkfile); + } +} +catch (const std::exception &e) { + std::cerr << "Program threw an exception: " << e.what() << "\n"; + throw; +} diff --git a/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp b/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp index 179fd8253..8e0b90f48 100644 --- a/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp +++ b/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp @@ -64,11 +64,11 @@ typedef Eigen::Array(fluid.numPhases() + 1, ADB::null())) + , residual_(std::vector(3, ADB::null())) { } @@ -82,7 +82,6 @@ typedef Eigen::Array& src, const std::vector& polymer_inflow) - // const bool if_polymer_actived) { V pvol(grid_.number_of_cells); @@ -100,6 +99,10 @@ typedef Eigen::Array(& x.saturation()[0], nc, np); for (int phase = 0; phase < np; ++phase) { state.saturation[phase] = ADB::constant(s_all.col(phase)); - // state.saturation[1] = ADB::constant(s_all.col(1)); } // Concentration - assert(not x.concentration().empty()); const V c = Eigen::Map(&x.concentration()[0], nc); @@ -256,29 +257,28 @@ typedef Eigen::Array(&polymer_inflow[0], nc); + const ADB src_polymer = polymerSource(kr ,src, polymer_inflow, state); ADB mc = computeMc(state); ADB mflux = computeMassFlux(0, trans, kr, state); residual_[2] = pvdt * (state.saturation[0] * state.concentration - old_state.saturation[0] * old_state.concentration) + ops_.div * state.concentration * mc * mflux - src_polymer; - } else { - residual_[2] = ADB::constant(V::Zero(nc)); - } + // } else { + // residual_[2] = ADB::constant(V::Zero(nc)); +// } + for (int i = 0; i < 3; ++i) + std::cout<<"residual_["<& polymer_inflow) const + polymerInjectedAmount(const std::vector& polymer_inflow) const { double amount = 0; for (int i = 0; i < int(polymer_inflow.size()); ++i) { @@ -292,7 +292,8 @@ typedef Eigen::Array& kr, - const std::vector& src) const + const std::vector& src + ) const { //extract the source to out and in source. std::vector outsrc; @@ -313,7 +314,6 @@ typedef Eigen::Array(& src[0], grid_.number_of_cells); const V outSrc = Eigen::Map(& outsrc[0], grid_.number_of_cells); const V inSrc = Eigen::Map(& insrc[0], grid_.number_of_cells); - // compute the out-fracflow. ADB f_out = computeFracFlow(phase, kr); // compute the in-fracflow. @@ -323,11 +323,44 @@ typedef Eigen::Array& kr, + const std::vector& src, + const std::vector& polymer_inflow_c, + const SolutionState& state) const + { + //extract the source to out and in source. + std::vector outsrc; + std::vector insrc; + std::vector::const_iterator it; + for (it = src.begin(); it != src.end(); ++it) { + if (*it < 0) { + outsrc.push_back(*it); + insrc.push_back(0.0); + } else if (*it > 0) { + insrc.push_back(*it); + outsrc.push_back(0.0); + } else { + outsrc.emplace_back(0); + insrc.emplace_back(0); + } + } + const V source = Eigen::Map(& src[0], grid_.number_of_cells); + const V outSrc = Eigen::Map(& outsrc[0], grid_.number_of_cells); + const V inSrc = Eigen::Map(& insrc[0], grid_.number_of_cells); + const V polyin = Eigen::Map(& polymer_inflow_c[0], grid_.number_of_cells); + // compute the out-fracflow. + ADB f_out = computeFracFlow(0, kr); + // compute the in-fracflow. + V f_in = V::Ones(grid_.number_of_cells); + return f_out * outSrc * state.concentration + f_in * inSrc * polyin ; + } ADB @@ -355,8 +388,8 @@ typedef Eigen::Array matr = mass_res.derivative()[0]; V dx(V::Zero(mass_res.size())); @@ -391,7 +424,6 @@ typedef Eigen::Array& kr, + const std::vector& src, + const std::vector& polymer_inflow_c, + const SolutionState& state) const; ADB computeMc(const SolutionState& state) const; @@ -100,7 +104,7 @@ namespace Opm { ADB transMult(const ADB& p) const; double - PolymerInjectedAmount(const std::vector& polymer_inflow) const; + polymerInjectedAmount(const std::vector& polymer_inflow) const; }; } // namespace Opm #endif// OPM_FULLYIMPLICITTWOPHASESOLVER_HEADER_INCLUDED diff --git a/opm/polymer/fullyimplicit/PolymerPropsAd.hpp b/opm/polymer/fullyimplicit/PolymerPropsAd.hpp index f03ccd34c..69936843d 100644 --- a/opm/polymer/fullyimplicit/PolymerPropsAd.hpp +++ b/opm/polymer/fullyimplicit/PolymerPropsAd.hpp @@ -65,7 +65,6 @@ namespace Opm { */ typedef AutoDiffBlock ADB; typedef ADB::V V; - PolymerPropsAd(const PolymerProperties& polymer_props); ~PolymerPropsAd();