OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2024-11-26 03:00:17 -06:00
Code Issues Packages Projects Releases Wiki Activity
a8f560b046
opm-simulators/tests/test_singlecellsolves.cpp
Atgeirr Flø Rasmussen cf9b7c39b9 Adapt to moved opm-grid headers.
2018-02-10 08:33:33 +01:00

278 lines
11 KiB
C++
Raw Blame History

/*
Copyright 2012 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif // HAVE_CONFIG_H
#include <opm/core/pressure/FlowBCManager.hpp>
#include <opm/grid/UnstructuredGrid.h>
#include <opm/grid/GridManager.hpp>
#include <opm/core/wells.h>
#include <opm/core/wells/WellsManager.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/core/simulator/initState.hpp>
#include <opm/core/simulator/SimulatorReport.hpp>
#include <opm/simulators/timestepping/SimulatorTimer.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/common/utility/parameters/ParameterGroup.hpp>
#include <opm/core/props/IncompPropertiesBasic.hpp>
#include <opm/core/props/IncompPropertiesFromDeck.hpp>
#include <opm/core/props/rock/RockCompressibility.hpp>
#include <opm/core/linalg/LinearSolverFactory.hpp>
#include <opm/polymer/PolymerState.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/polymer/IncompTpfaPolymer.hpp>
#include <opm/polymer/TransportSolverTwophasePolymer.hpp>
#include <opm/polymer/PolymerProperties.hpp>
#include <boost/scoped_ptr.hpp>
#include <algorithm>
#include <iostream>
#include <vector>
#include <numeric>
// ----------------- Main program -----------------
int
main(int argc, char** argv)
try
{
using namespace Opm;
// std::cout << "\n================ Test program for single-cell solves with polymer ===============\n\n";
ParameterGroup param(argc, argv, false);
param.disableOutput();
// std::cout << "--------------- Reading parameters ---------------" << std::endl;
boost::scoped_ptr<GridManager> grid;
boost::scoped_ptr<IncompPropertiesInterface> props;
std::unique_ptr<PolymerState> state;
Opm::PolymerProperties poly_props;
// bool check_well_controls = false;
// int max_well_control_iterations = 0;
// -------- Initialising section ----------
// Grid init.
grid.reset(new GridManager(2, 1, 1, 1.0, 1.0, 1.0));
// Rock and fluid init.
{
const UnstructuredGrid& ug_grid = *(grid->c_grid());
state.reset( new PolymerState( UgGridHelpers::numCells( ug_grid ) , UgGridHelpers::numFaces( ug_grid ), 2));
props.reset(new IncompPropertiesBasic(param, ug_grid.dimensions, UgGridHelpers::numCells( ug_grid )));
// Init state variables (saturation and pressure).
initStateBasic(*grid->c_grid(), *props, param, 0.0, *state);
// Init Polymer state
if (param.has("poly_init")) {
double poly_init = param.getDefault("poly_init", 0.0);
for (int cell = 0; cell < grid->c_grid()->number_of_cells; ++cell) {
double smin[2], smax[2];
props->satRange(1, &cell, smin, smax);
auto& saturation = state->saturation();
auto& concentration = state->getCellData( state->CONCENTRATION );
auto& max_concentration = state->getCellData( state->CMAX );
props->satRange(1, &cell, smin, smax);
if (saturation[2*cell] > 0.5*(smin[0] + smax[0])) {
concentration[cell] = poly_init;
max_concentration[cell] = poly_init;
} else {
saturation[2*cell + 0] = 0.;
saturation[2*cell + 1] = 1.;
concentration[cell] = 0.;
max_concentration[cell] = 0.;
}
}
}
}
// Init polymer properties.
// Setting defaults to provide a simple example case.
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.0); // Note that we default to no dps here!
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<int>("ads_index", Opm::PolymerProperties::NoDesorption);
std::vector<double> c_vals_visc(2, -1e100);
c_vals_visc[0] = 0.0;
c_vals_visc[1] = 7.0;
std::vector<double> 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<double> 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<double> ads_vals(3, -1e100);
ads_vals[0] = 0.0;
ads_vals[1] = 0.0015;
ads_vals[2] = 0.0025;
// ads_vals[1] = 0.0;
// ads_vals[2] = 0.0;
std::vector<double> water_vel_vals(2, -1e100);
water_vel_vals[0] = 0.0;
water_vel_vals[1] = 10.0;
std::vector<double> shear_vrf_vals(2, -1e100);
shear_vrf_vals[0] = 1.0;
shear_vrf_vals[1] = 1.0;
poly_props.set(c_max, mix_param, rock_density, dead_pore_vol, res_factor, c_max_ads,
static_cast<Opm::PolymerProperties::AdsorptionBehaviour>(ads_index),
c_vals_visc, visc_mult_vals, c_vals_ads, ads_vals, water_vel_vals, shear_vrf_vals);
// Initialising src
int num_cells = grid->c_grid()->number_of_cells;
std::vector<double> src(num_cells, 0.0);
// Compute pore volumes, in order to enable specifying injection rate
// terms of total pore volume.
std::vector<double> porevol;
computePorevolume(*grid->c_grid(), props->porosity(), porevol);
const double default_injection = 1.0;
const double flow_per_sec = param.getDefault<double>("injected_porevolumes_per_sec", default_injection)
*porevol[0];
src[0] = flow_per_sec;
src[num_cells - 1] = -flow_per_sec;
// Boundary conditions.
FlowBCManager bcs;
// Linear solver.
LinearSolverFactory linsolver(param);
// Reordering solver.
const double nl_tolerance = param.getDefault("nl_tolerance", 1e-9);
const int nl_maxiter = param.getDefault("nl_maxiter", 30);
Opm::TransportSolverTwophasePolymer::SingleCellMethod method;
std::string method_string = param.getDefault("single_cell_method", std::string("Bracketing"));
if (method_string == "Bracketing") {
method = Opm::TransportSolverTwophasePolymer::Bracketing;
} else if (method_string == "Newton") {
method = Opm::TransportSolverTwophasePolymer::Newton;
} else if (method_string == "Gradient") {
method = Opm::TransportSolverTwophasePolymer::Gradient;
} else if (method_string == "NewtonSimpleSC") {
method = Opm::TransportSolverTwophasePolymer::NewtonSimpleSC;
} else if (method_string == "NewtonSimpleC") {
method = Opm::TransportSolverTwophasePolymer::NewtonSimpleC;
} else {
OPM_THROW(std::runtime_error, "Unknown method: " << method_string);
}
Opm::TransportSolverTwophasePolymer reorder_model(*grid->c_grid(), *props, poly_props,
method, nl_tolerance, nl_maxiter);
// Warn if any parameters are unused.
// if (param.anyUnused()) {
// std::cout << "-------------------- Unused parameters: --------------------\n";
// param.displayUsage();
// std::cout << "----------------------------------------------------------------" << std::endl;
// }
// Write parameters to file for later reference.
param.writeParam("test_singlecellsolves.param");
// Setting up a number of input (s, c) pairs and solving.
// HACK warning: we manipulate the source term,
// but the compressibility term in the solver
// assumes that all inflow is water inflow. Therefore
// one must zero the compressibility term in
// TransportSolverTwophasePolymer line 365 before compiling this program.
// (To fix this we should add proper all-phase src terms.)
std::vector<double> transport_src = src;
const double dt = param.getDefault("dt", 1.0);
const int num_sats = 501;
const int num_concs = 501;
// Find the face between cell 0 and 1...
const UnstructuredGrid& ug = *grid->c_grid();
int face01 = -1;
for (int f = 0; f < ug.number_of_faces; ++f) {
if (ug.face_cells[2*f] == 0 && ug.face_cells[2*f+1] == 1) {
face01 = f;
break;
}
}
if (face01 == -1) {
OPM_THROW(std::runtime_error, "Could not find face adjacent to cells [0 1]");
}
state->faceflux()[face01] = src[0];
for (int sats = 0; sats < num_sats; ++sats) {
const double s = double(sats)/double(num_sats - 1);
const double ff = s; // Simplified a lot...
for (int conc = 0; conc < num_concs; ++conc) {
const double c = poly_props.cMax()*double(conc)/double(num_concs - 1);
std::vector<double> polymer_inflow_c(num_cells, c);
// std::cout << "(s, c) = (" << s << ", " << c << ")\n";
transport_src[0] = src[0]*ff;
// Resetting the state for next run.
auto& saturation = state->saturation();
auto& concentration = state->getCellData( state->CONCENTRATION );
auto& max_concentration = state->getCellData( state->CMAX );
saturation[0] = 0.0;
saturation[1] = 0.0;
concentration[0] = 0.0;
concentration[1] = 0.0;
max_concentration[0] = 0.0;
max_concentration[1] = 0.0;
reorder_model.solve(&state->faceflux()[0],
&porevol[0],
&transport_src[0],
&polymer_inflow_c[0],
dt,
saturation,
concentration,
max_concentration);
#ifdef PROFILING
// Extract residual counts.
typedef std::list<Opm::TransportSolverTwophasePolymer::Newton_Iter> ListRes;
const ListRes& res_counts = reorder_model.res_counts;
double counts[2] = { 0, 0 };
for (ListRes::const_iterator it = res_counts.begin(); it != res_counts.end(); ++it) {
if (it->cell == 0) {
++counts[it->res_s];
}
}
// std::cout << "c residual count: " << counts[0] << '\n';
// std::cout << "s residual count: " << counts[1] << '\n';
std::cout << counts[0] << ' ' << counts[1] << ' ' << s << ' ' << c << '\n';
#endif
}
}
}
catch (const std::exception &e) {
std::cerr << "Program threw an exception: " << e.what() << "\n";
throw;
}
Reference in New Issue View Git Blame Copy Permalink