diff --git a/examples/polymer_reorder.cpp b/examples/polymer_reorder.cpp index 55ee126b6..f857605fd 100644 --- a/examples/polymer_reorder.cpp +++ b/examples/polymer_reorder.cpp @@ -237,6 +237,10 @@ public: return smax_[2*c + 0]; } + double cMax() const + { + return polyprops_.cMax(); + } template NewtonPolymerTransportModel; +typedef Opm::SinglePointUpwindTwoPhasePolymer FluxModel; using namespace Opm::ImplicitTransportDefault; @@ -271,7 +275,7 @@ public: } }; -typedef Opm::ImplicitTransportc_grid()->number_of_cells; ++cell) { - if (state.saturation()[2*cell] > 0) { + double smin[2], smax[2]; + props->satRange(1, &cell, smin, smax); + if (state.saturation()[2*cell] > 0.5*(smin[0] + smax[0])) { state.concentration()[cell] = poly_init; state.maxconcentration()[cell] = poly_init; } else { + state.saturation()[2*cell + 0] = 0.; + state.saturation()[2*cell + 1] = 1.; state.concentration()[cell] = 0.; - state.maxconcentration()[cell] = poly_init; + state.maxconcentration()[cell] = 0.; } } } @@ -435,10 +443,10 @@ main(int argc, char** argv) std::vector ads_vals(3, -1e100); ads_vals[0] = 0.0; // polyprop.ads_vals[1] = param.getDefault("c_max_ads", 0.0025); - // ads_vals[1] = 0.0015; - // ads_vals[2] = 0.0025; - ads_vals[1] = 0.0; - ads_vals[2] = 0.0; + ads_vals[1] = 0.0015; + ads_vals[2] = 0.0025; + // ads_vals[1] = 0.0; + // ads_vals[2] = 0.0; polyprop.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); @@ -552,17 +560,17 @@ main(int argc, char** argv) reorder_model.initGravity(grav); } // Non-reordering solver. - NewtonPolymerTransportModel model(fluid, *grid->c_grid(), porevol, grav, guess_old_solution); + FluxModel fmodel(fluid, *grid->c_grid(), porevol, grav, guess_old_solution); if (use_gravity) { - model.initGravityTrans(*grid->c_grid(), psolver.getHalfTrans()); + fmodel.initGravityTrans(*grid->c_grid(), psolver.getHalfTrans()); } - TransportSolver tsolver(model); + TransportSolver tsolver(fmodel); // Column-based gravity segregation solver. std::vector > columns; if (use_column_solver) { Opm::extractColumn(*grid->c_grid(), columns); } - Opm::GravityColumnSolverPolymer colsolver(model, *grid->c_grid(), nl_tolerance, nl_maxiter); + Opm::GravityColumnSolverPolymer colsolver(fmodel, fluid, *grid->c_grid(), nl_tolerance, nl_maxiter); // // // Not implemented for polymer. // // Control init. diff --git a/opm/polymer/GravityColumnSolverPolymer.hpp b/opm/polymer/GravityColumnSolverPolymer.hpp index e9f6f86f8..922c0a46b 100644 --- a/opm/polymer/GravityColumnSolverPolymer.hpp +++ b/opm/polymer/GravityColumnSolverPolymer.hpp @@ -29,16 +29,17 @@ namespace Opm /// Class for doing gravity segregation (only), /// on a vertical column of cells. - template + template class GravityColumnSolverPolymer { public: /// Note: the model will be changed since it stores computed /// quantities in itself, such as mobilities. - GravityColumnSolverPolymer(Model& model, - const UnstructuredGrid& grid, - const double tol, - const int maxit); + GravityColumnSolverPolymer(FluxModel& fmodel, + const Model& model, + const UnstructuredGrid& grid, + const double tol, + const int maxit); /// \param[in] columns for each column (with logical cartesian indices as key), /// contains the cells on which to solve the segregation @@ -59,7 +60,8 @@ namespace Opm std::vector& cmax, std::vector& sol_vec ); - Model& model_; + FluxModel& fmodel_; + const Model& model_; const UnstructuredGrid& grid_; const double tol_; const int maxit_; diff --git a/opm/polymer/GravityColumnSolverPolymer_impl.hpp b/opm/polymer/GravityColumnSolverPolymer_impl.hpp index efb67f2f7..c47275a17 100644 --- a/opm/polymer/GravityColumnSolverPolymer_impl.hpp +++ b/opm/polymer/GravityColumnSolverPolymer_impl.hpp @@ -41,15 +41,24 @@ #include #include +// namespace { + +// int sgn(double val) { +// return (0.0 < val) - (val < 0); +// } + +// } + namespace Opm { - template - GravityColumnSolverPolymer::GravityColumnSolverPolymer(Model& model, - const UnstructuredGrid& grid, - const double tol, - const int maxit) - : model_(model), grid_(grid), tol_(tol), maxit_(maxit) + template + GravityColumnSolverPolymer::GravityColumnSolverPolymer(FluxModel& fmodel, + const Model& model, + const UnstructuredGrid& grid, + const double tol, + const int maxit) + : fmodel_(fmodel), model_(model), grid_(grid), tol_(tol), maxit_(maxit) { } @@ -117,8 +126,8 @@ namespace Opm /// problem. For each column, its cells must be in a single /// vertical column, connected and ordered /// (direction doesn't matter). - template - void GravityColumnSolverPolymer::solve(const std::vector >& columns, + template + void GravityColumnSolverPolymer::solve(const std::vector >& columns, const double dt, std::vector& s, std::vector& c, @@ -129,20 +138,83 @@ namespace Opm StateWithZeroFlux state(s, c, cmax); // This holds s, c and cmax by reference. JacSys sys(2*grid_.number_of_cells); std::vector increment(2*grid_.number_of_cells, 0.0); - model_.initStep(state, grid_, sys); + fmodel_.initStep(state, grid_, sys); int iter = 0; double max_delta = 1e100; + const double cmax_cell = 2.0*model_.cMax(); + const double tol_c_cell = 1e-2*cmax_cell; while (iter < maxit_) { - model_.initIteration(state, grid_, sys); + fmodel_.initIteration(state, grid_, sys); int size = columns.size(); for(int i = 0; i < size; ++i) { solveSingleColumn(columns[i], dt, s, c, cmax, increment); } for (int cell = 0; cell < grid_.number_of_cells; ++cell) { - sys.vector().writableSolution()[2*cell + 0] += increment[2*cell + 0]; - sys.vector().writableSolution()[2*cell + 1] += increment[2*cell + 1]; - } + double& s_cell = sys.vector().writableSolution()[2*cell + 0]; + double& c_cell = sys.vector().writableSolution()[2*cell + 1]; + s_cell += increment[2*cell + 0]; + c_cell += increment[2*cell + 1]; + if (s_cell < 0.) { + double& incr = increment[2*cell + 0]; + s_cell -= incr; + if (std::fabs(incr) < 1e-2) { + incr = -s_cell; + s_cell = 0.; + } else { + incr = -s_cell/2.0; + s_cell = s_cell/2.0; + } + } + if (s_cell > 1.) { + double& incr = increment[2*cell + 0]; + s_cell -= incr; + if (std::fabs(incr) < 1e-2) { + incr = 1. - s_cell; + s_cell = 1.; + } else { + incr = (1 - s_cell)/2.0; + s_cell = (1 + s_cell)/2.0; + } + } + if (c_cell < 0.) { + double& incr = increment[2*cell + 1]; + c_cell -= incr; + if (std::fabs(incr) < tol_c_cell) { + incr = -c_cell; + c_cell = 0.; + } else { + incr = -c_cell/2.0; + c_cell = c_cell/2.0; + } + } + if (c_cell > cmax_cell) { + double& incr = increment[2*cell + 1]; + c_cell -= incr; + if (std::fabs(incr) < tol_c_cell) { + incr = cmax_cell - c_cell; + c_cell = cmax_cell; + } else { + incr = (cmax_cell - c_cell)/2.0; + c_cell = (cmax_cell + c_cell)/2.0; + } + } + + // if (s_cell < 0.) { + // increment[2*cell + 0] = increment[2*cell + 0] - s_cell; + // s_cell = 0.; + // } else if (s_cell > 1.) { + // increment[2*cell + 0] = increment[2*cell + 0] - s_cell + 1.; + // s_cell = 1.; + // } + // if (c_cell < 0) { + // increment[2*cell + 1] = increment[2*cell + 1] - c_cell; + // c_cell = 0.; + // } else if (c_cell > cmax_cell) { + // increment[2*cell + 1] = increment[2*cell + 1] - c_cell + cmax_cell; + // c_cell = cmax_cell; + // } + } const double maxelem = *std::max_element(increment.begin(), increment.end()); const double minelem = *std::min_element(increment.begin(), increment.end()); max_delta = std::max(maxelem, -minelem); @@ -156,10 +228,10 @@ namespace Opm THROW("Failed to converge!"); } // Finalize. - // model_.finishIteration(); // + // fmodel_.finishIteration(); // // finishStep() writes to state, which holds s by reference. // This will update the entire grid's state... cmax is updated here. - model_.finishStep(grid_, sys.vector().solution(), state); + fmodel_.finishStep(grid_, sys.vector().solution(), state); } @@ -168,8 +240,8 @@ namespace Opm /// \param[in] column_cells the cells on which to solve the segregation /// problem. Must be in a single vertical column, /// and ordered (direction doesn't matter). - template - void GravityColumnSolverPolymer::solveSingleColumn(const std::vector& column_cells, + template + void GravityColumnSolverPolymer::solveSingleColumn(const std::vector& column_cells, const double dt, std::vector& s, std::vector& c, @@ -215,7 +287,7 @@ namespace Opm F.assign(2, 0.); dFd1.assign(4, 0.); dFd2.assign(4, 0.); - model_.fluxConnection(state, grid_, dt, cell, face, &F[0], &dFd1[0], &dFd2[0]); + fmodel_.fluxConnection(state, grid_, dt, cell, face, &F[0], &dFd1[0], &dFd2[0]); if (c1 == prev_cell || c2 == prev_cell) { hm[bmc(2*ci + 0, 2*(ci - 1) + 0)] += dFd2[0]; hm[bmc(2*ci + 0, 2*(ci - 1) + 1)] += dFd2[1]; @@ -239,7 +311,7 @@ namespace Opm } F.assign(2, 0.); dF.assign(4, 0.); - model_.accumulation(grid_, cell, &F[0], &dF[0]); + fmodel_.accumulation(grid_, cell, &F[0], &dF[0]); hm[bmc(2*ci + 0, 2*ci + 0)] += dF[0]; hm[bmc(2*ci + 0, 2*ci + 1)] += dF[1]; hm[bmc(2*ci + 1, 2*ci + 0)] += dF[2];