OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge from upstream.

Browse Source
This commit is contained in:
Halvor M. Nilsen 2012-06-12 12:47:45 +02:00
commit 97415da35d
6 changed files with 77 additions and 305 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

28
generate_doc_figures.py
View File

@ -90,6 +90,34 @@ for case in cases:
WriteImage(join(figure_path, "tutorial3-"+case+".png")) WriteImage(join(figure_path, "tutorial3-"+case+".png"))
Hide(grid) Hide(grid)
# tutorial 4
call(join(tutorial_path, "tutorial4"))
for case in range(0,20):
data_file_name = join(tutorial_data_path, "tutorial4-"+"%(case)03d"%{"case": case}+".vtu")
collected_garbage_file.append(data_file_name)
cases = ["000", "005", "010", "015", "019"]
for case in cases:
data_file_name = join(tutorial_data_path, "tutorial4-"+case+".vtu")
grid = XMLUnstructuredGridReader(FileName = data_file_name)
grid.UpdatePipeline()
Show(grid)
dp = GetDisplayProperties(grid)
dp.Representation = 'Surface'
dp.ColorArrayName = 'saturation'
sat = grid.CellData.GetArray(1)
sat_lookuptable = GetLookupTableForArray( "saturation", 1, RGBPoints=[0, 1, 0, 0, 1, 0, 0, 1])
dp.LookupTable = sat_lookuptable
view.Background = [1, 1, 1]
camera = GetActiveCamera()
camera.SetPosition(100, 100, 550)
camera.SetViewUp(0, 1, 0)
camera.SetViewAngle(30)
camera.SetFocalPoint(100, 100, 5)
Render()
WriteImage(join(figure_path, "tutorial4-"+case+".png"))
Hide(grid)
# remove temporary files # remove temporary files
for f in collected_garbage_file: for f in collected_garbage_file:
remove(f) remove(f)

16
opm/core/grid.h
View File

@ -250,18 +250,20 @@ struct UnstructuredGrid *
create_grid_empty(void); create_grid_empty(void);
/** /**
Allocate and initialise an UnstructuredGrid where pointers are set to Allocate and initialise an UnstructuredGrid where pointers are set
location with correct size. to location with correct size.
\param[in] ndims Number of physical dimensions. \param[in] ndims Number of physical dimensions.
\param[in] ncells Number of cells. \param[in] ncells Number of cells.
\param[in] nfaces Number of faces. \param[in] nfaces Number of faces.
\param[in] nfacenodes Size of mapping from faces to nodes. \param[in] nfacenodes Size of mapping from faces to nodes.
\param[in] ncellfaces Size of mapping from cells to faces.(i.e Number of halffaces) \param[in] ncellfaces Size of mapping from cells to faces.
\param[in] nnodes Number of Nodes (i.e., the number of `half-faces')
\return Fully formed UnstructuredGrid with all fields allocated, but not filled with \param[in] nnodes Number of Nodes.
usefull numbers.
<code>NULL</code> in case of allocation failure. \return Fully formed UnstructuredGrid with all fields except
<code>global_cell</code> allocated, but not filled with meaningful
values. <code>NULL</code> in case of allocation failure.
*/ */
struct UnstructuredGrid * struct UnstructuredGrid *
allocate_grid(size_t ndims , allocate_grid(size_t ndims ,

65
opm/core/grid/cart_grid.c
View File

@ -170,68 +170,13 @@ allocate_cart_grid(size_t ndims ,
size_t nfaces, size_t nfaces,
size_t nnodes) size_t nnodes)
{ {
size_t nel; size_t nfacenodes, ncellfaces;
struct UnstructuredGrid *G;
G = create_grid_empty(); nfacenodes = nfaces * (2 * (ndims - 1));
ncellfaces = ncells * (2 * ndims);
if (G != NULL) { return allocate_grid(ndims, ncells, nfaces,
/* Node fields ---------------------------------------- */ nfacenodes, ncellfaces, nnodes);
nel = nnodes * ndims;
G->node_coordinates = malloc(nel * sizeof *G->node_coordinates);
/* Face fields ---------------------------------------- */
nel = nfaces * (2 * (ndims - 1));
G->face_nodes = malloc(nel * sizeof *G->face_nodes);
nel = nfaces + 1;
G->face_nodepos = malloc(nel * sizeof *G->face_nodepos);
nel = 2 * nfaces;
G->face_cells = malloc(nel * sizeof *G->face_cells);
nel = nfaces * ndims;
G->face_centroids = malloc(nel * sizeof *G->face_centroids);
nel = nfaces * ndims;
G->face_normals = malloc(nel * sizeof *G->face_normals);
nel = nfaces * 1;
G->face_areas = malloc(nel * sizeof *G->face_areas);
/* Cell fields ---------------------------------------- */
nel = ncells * (2 * ndims);
G->cell_faces = malloc(nel * sizeof *G->cell_faces);
nel = ncells + 1;
G->cell_facepos = malloc(nel * sizeof *G->cell_facepos);
nel = ncells * ndims;
G->cell_centroids = malloc(nel * sizeof *G->cell_centroids);
nel = ncells * 1;
G->cell_volumes = malloc(nel * sizeof *G->cell_volumes);
if ((G->node_coordinates == NULL) ||
(G->face_nodes == NULL) ||
(G->face_nodepos == NULL) ||
(G->face_cells == NULL) ||
(G->face_centroids == NULL) ||
(G->face_normals == NULL) ||
(G->face_areas == NULL) ||
(G->cell_faces == NULL) ||
(G->cell_facepos == NULL) ||
(G->cell_centroids == NULL) ||
(G->cell_volumes == NULL) )
{
destroy_grid(G);
G = NULL;
}
}
return G;
} }

21
opm/core/linalg/LinearSolverFactory.cpp
View File

@ -55,20 +55,29 @@ namespace Opm
LinearSolverFactory::LinearSolverFactory(const parameter::ParameterGroup& param) LinearSolverFactory::LinearSolverFactory(const parameter::ParameterGroup& param)
{ {
const std::string ls = param.getDefault<std::string>("linsolver", "umfpack"); const std::string ls =
param.getDefault<std::string>("linsolver", "umfpack");
if (ls == "umfpack") { if (ls == "umfpack") {
#if HAVE_SUITESPARSE_UMFPACK_H #if HAVE_SUITESPARSE_UMFPACK_H
solver_.reset(new LinearSolverUmfpack); solver_.reset(new LinearSolverUmfpack);
#else
THROW("Linear solver " << ls <<" is not available.");
#endif #endif
} else if (ls == "istl") { }
else if (ls == "istl") {
#if HAVE_DUNE_ISTL #if HAVE_DUNE_ISTL
solver_.reset(new LinearSolverIstl(param)); solver_.reset(new LinearSolverIstl(param));
#else
THROW("Linear solver " << ls <<" is not available.");
#endif #endif
} }
else {
THROW("Linear solver " << ls << " is unknown.");
}
if (! solver_) {
THROW("Linear solver " << ls << " is not enabled in "
"this configuration.");
}
} }

1
opm/core/wells/WellsGroup.cpp
View File

@ -941,7 +941,6 @@ namespace Opm
double WellNode::productionGuideRate(bool only_group) double WellNode::productionGuideRate(bool only_group)
{ {
if (!only_group || prodSpec().control_mode_ == ProductionSpecification::GRUP) { if (!only_group || prodSpec().control_mode_ == ProductionSpecification::GRUP) {
std::cout << prodSpec().guide_rate_ << std::endl;
return prodSpec().guide_rate_; return prodSpec().guide_rate_;
} }
return 0.0; return 0.0;

251
tutorials/tutorial4.cpp
View File

@ -35,14 +35,7 @@
#include <opm/core/pressure/FlowBCManager.hpp> #include <opm/core/pressure/FlowBCManager.hpp>
#include <opm/core/fluid/IncompPropertiesBasic.hpp> #include <opm/core/fluid/IncompPropertiesBasic.hpp>
#include <opm/core/transport/transport_source.h> #include <opm/core/transport/reorder/TransportModelTwophase.hpp>
#include <opm/core/transport/CSRMatrixUmfpackSolver.hpp>
#include <opm/core/transport/NormSupport.hpp>
#include <opm/core/transport/ImplicitAssembly.hpp>
#include <opm/core/transport/ImplicitTransport.hpp>
#include <opm/core/transport/JacobianSystem.hpp>
#include <opm/core/transport/CSRMatrixBlockAssembler.hpp>
#include <opm/core/transport/SinglePointUpwindTwoPhase.hpp>
#include <opm/core/simulator/TwophaseState.hpp> #include <opm/core/simulator/TwophaseState.hpp>
@ -52,113 +45,6 @@
#include <opm/core/wells/WellCollection.hpp> #include <opm/core/wells/WellCollection.hpp>
/// \page tutorial4 Well controls /// \page tutorial4 Well controls
///
/// \page tutorial4
/// \section commentedcode4 Commented code:
/// \page tutorial4
/// \details
/// We define a class which computes mobility, capillary pressure and
/// the minimum and maximum saturation value for each cell.
/// \code
class TwophaseFluid
{
public:
/// \endcode
/// \page tutorial4
/// \details Constructor operator. Takes in the fluid properties defined
/// \c props
/// \code
TwophaseFluid(const Opm::IncompPropertiesInterface& props);
/// \endcode
/// \page tutorial4
/// \details Density for each phase.
/// \code
double density(int phase) const;
/// \endcode
/// \page tutorial4
/// \details Computes the mobility and its derivative with respect to saturation
/// for a given cell \c c and saturation \c sat. The template classes \c Sat,
/// \c Mob, \c DMob are typically arrays. By using templates, we do not have to
/// investigate how these array objects are implemented
/// (as long as they have an \c operator[] method).
/// \code
template <class Sat, class Mob, class DMob>
void mobility(int c, const Sat& s, Mob& mob, DMob& dmob) const;
/// \endcode
/// \page tutorial4
/// \details Computes the capillary pressure and its derivative with respect
/// to saturation
/// for a given cell \c c and saturation \c sat.
/// \code
template <class Sat, class Pcap, class DPcap>
void pc(int c, const Sat& s, Pcap& pcap, DPcap& dpcap) const;
/// \endcode
/// \page tutorial4
/// \details Returns the minimum permitted saturation.
/// \code
double s_min(int c) const;
/// \endcode
/// \page tutorial4
/// \details Returns the maximum permitted saturation
/// \code
double s_max(int c) const;
/// \endcode
/// \page tutorial4
/// \details Private variables
/// \code
private:
const Opm::IncompPropertiesInterface& props_;
std::vector<double> smin_;
std::vector<double> smax_;
};
/// \endcode
/// \page tutorial4
/// \details We set up the transport model.
/// \code
typedef Opm::SinglePointUpwindTwoPhase<TwophaseFluid> TransportModel;
/// \endcode
/// \page tutorial4
/// \details
/// The transport equation is nonlinear. We use an implicit transport solver
/// which implements a Newton-Raphson solver.
/// We define the format of the objects
/// which will be used by the solver.
/// \code
using namespace Opm::ImplicitTransportDefault;
typedef NewtonVectorCollection< ::std::vector<double> > NVecColl;
typedef JacobianSystem< struct CSRMatrix, NVecColl > JacSys;
template <class Vector>
class MaxNorm {
public:
static double
norm(const Vector& v) {
return AccumulationNorm <Vector, MaxAbs>::norm(v);
}
};
/// \endcode
/// \page tutorial4
/// \details
/// We set up the solver.
/// \code
typedef Opm::ImplicitTransport<TransportModel,
JacSys ,
MaxNorm ,
VectorNegater ,
VectorZero ,
MatrixZero ,
VectorAssign > TransportSolver;
/// \endcode
/// \page tutorial4 /// \page tutorial4
/// \details /// \details
@ -179,8 +65,9 @@ int main ()
double dy = 10.; double dy = 10.;
double dz = 10.; double dz = 10.;
using namespace Opm; using namespace Opm;
GridManager grid(nx, ny, nz, dx, dy, dz); GridManager grid_manager(nx, ny, nz, dx, dy, dz);
int num_cells = grid.c_grid()->number_of_cells; const UnstructuredGrid& grid = *grid_manager.c_grid();
int num_cells = grid.number_of_cells;
/// \endcode /// \endcode
/// \page tutorial4 /// \page tutorial4
@ -214,8 +101,7 @@ int main ()
/// description and set this function to be linear. For more realistic fluid, the /// description and set this function to be linear. For more realistic fluid, the
/// saturation function is given by the data. /// saturation function is given by the data.
/// \code /// \code
SaturationPropsBasic::RelPermFunc rel_perm_func; SaturationPropsBasic::RelPermFunc rel_perm_func = SaturationPropsBasic::Linear;
rel_perm_func = SaturationPropsBasic::Linear;
/// \endcode /// \endcode
/// \page tutorial4 /// \page tutorial4
@ -224,7 +110,6 @@ int main ()
/// \code /// \code
IncompPropertiesBasic props(num_phases, rel_perm_func, rho, mu, IncompPropertiesBasic props(num_phases, rel_perm_func, rho, mu,
porosity, k, dim, num_cells); porosity, k, dim, num_cells);
TwophaseFluid fluid(props);
/// \endcode /// \endcode
@ -237,15 +122,6 @@ int main ()
/// \page tutorial4
/// \details We set up the source term
/// \code
std::vector<double> src(num_cells, 0.0);
src[0] = 1.;
src[num_cells-1] = -1.;
/// \endcode
/// \page tutorial4 /// \page tutorial4
/// \details We set up necessary information for the wells /// \details We set up necessary information for the wells
/// \code /// \code
@ -259,34 +135,27 @@ int main ()
/// \endcode /// \endcode
/// \page tutorial4 /// \page tutorial4
/// \details We set up the source term for the transport solver. /// \details We set up the source term. Positive numbers indicate that the cell is a source,
/// while negative numbers indicate a sink.
/// \code /// \code
TransportSource* tsrc = create_transport_source(2, 2); std::vector<double> src(num_cells, 0.0);
double ssrc[] = { 1.0, 0.0 }; src[0] = 1.;
double ssink[] = { 0.0, 1.0 }; src[num_cells-1] = -1.;
double zdummy[] = { 0.0, 0.0 };
for (int cell = 0; cell < num_cells; ++cell) {
if (src[cell] > 0.0) {
append_transport_source(cell, 2, 0, src[cell], ssrc, zdummy, tsrc);
} else if (src[cell] < 0.0) {
append_transport_source(cell, 2, 0, src[cell], ssink, zdummy, tsrc);
}
}
/// \endcode /// \endcode
/// \page tutorial4 /// \page tutorial4
/// \details We compute the pore volume /// \details We compute the pore volume
/// \code /// \code
std::vector<double> porevol; std::vector<double> porevol;
Opm::computePorevolume(*grid.c_grid(), props.porosity(), porevol); Opm::computePorevolume(grid, props.porosity(), porevol);
/// \endcode /// \endcode
/// \page tutorial4 /// \page tutorial4
/// \details We set up the transport solver. /// \details Set up the transport solver. This is a reordering implicit Euler transport solver.
/// \code /// \code
const bool guess_old_solution = true; const double tolerance = 1e-9;
TransportModel model (fluid, *grid.c_grid(), porevol, grav, guess_old_solution); const int max_iterations = 30;
TransportSolver tsolver(model); Opm::TransportModelTwophase transport_solver(grid, props, tolerance, max_iterations);
/// \endcode /// \endcode
/// \page tutorial4 /// \page tutorial4
@ -296,14 +165,6 @@ int main ()
int num_time_steps = 20; int num_time_steps = 20;
/// \endcode /// \endcode
/// \page tutorial4
/// \details Control paramaters for the implicit solver.
/// \code
ImplicitTransportDetails::NRReport rpt;
ImplicitTransportDetails::NRControl ctrl;
/// \endcode
/// \page tutorial4 /// \page tutorial4
/// \details We define a vector which contains all cell indexes. We use this /// \details We define a vector which contains all cell indexes. We use this
@ -321,21 +182,13 @@ int main ()
FlowBCManager bcs; FlowBCManager bcs;
/// \endcode /// \endcode
/// \page tutorial4
/// \details
/// Linear solver init.
/// \code
using ImplicitTransportLinAlgSupport::CSRMatrixUmfpackSolver;
CSRMatrixUmfpackSolver linsolve;
/// \endcode
/// \page tutorial4 /// \page tutorial4
/// \details /// \details
/// We set up a two-phase state object, and /// We set up a two-phase state object, and
/// initialise water saturation to minimum everywhere. /// initialise water saturation to minimum everywhere.
/// \code /// \code
TwophaseState state; TwophaseState state;
state.init(*grid.c_grid(), 2); state.init(grid, 2);
state.setFirstSat(allcells, props, TwophaseState::MinSat); state.setFirstSat(allcells, props, TwophaseState::MinSat);
/// \endcode /// \endcode
@ -450,7 +303,7 @@ int main ()
/// last argument. /// last argument.
/// \code /// \code
LinearSolverUmfpack linsolver; LinearSolverUmfpack linsolver;
IncompTpfa psolver(*grid.c_grid(), props.permeability(), grav, linsolver, wells); IncompTpfa psolver(grid, props.permeability(), grav, linsolver, wells);
/// \endcode /// \endcode
@ -469,7 +322,7 @@ int main ()
/// \page tutorial4 /// \page tutorial4
/// \details In order to use the well controls, we need to generate the WDP for each well. /// \details In order to use the well controls, we need to generate the WDP for each well.
/// \code /// \code
Opm::computeWDP(*wells, *grid.c_grid(), state.saturation(), props.density(), gravity, true, wdp); Opm::computeWDP(*wells, grid, state.saturation(), props.density(), gravity, true, wdp);
/// \endcode /// \endcode
/// \page tutorial4 /// \page tutorial4
@ -511,7 +364,7 @@ int main ()
/// \page tutorial4 /// \page tutorial4
/// \details Transport solver /// \details Transport solver
/// \code /// \code
tsolver.solve(*grid.c_grid(), tsrc, dt, ctrl, state, linsolve, rpt); transport_solver.solve(&state.faceflux()[0], &porevol[0], &src[0], dt, state.saturation());
/// \endcode /// \endcode
/// \page tutorial4 /// \page tutorial4
@ -523,75 +376,11 @@ int main ()
Opm::DataMap dm; Opm::DataMap dm;
dm["saturation"] = &state.saturation(); dm["saturation"] = &state.saturation();
dm["pressure"] = &state.pressure(); dm["pressure"] = &state.pressure();
Opm::writeVtkData(*grid.c_grid(), dm, vtkfile); Opm::writeVtkData(grid, dm, vtkfile);
} }
} }
/// \endcode /// \endcode
/// \page tutorial4
/// \details Implementation of the TwophaseFluid class.
/// \code
TwophaseFluid::TwophaseFluid(const Opm::IncompPropertiesInterface& props)
: props_(props),
smin_(props.numCells()*props.numPhases()),
smax_(props.numCells()*props.numPhases())
{
const int num_cells = props.numCells();
std::vector<int> cells(num_cells);
for (int c = 0; c < num_cells; ++c) {
cells[c] = c;
}
props.satRange(num_cells, &cells[0], &smin_[0], &smax_[0]);
}
double TwophaseFluid::density(int phase) const
{
return props_.density()[phase];
}
template <class Sat,
class Mob,
class DMob>
void TwophaseFluid::mobility(int c, const Sat& s, Mob& mob, DMob& dmob) const
{
props_.relperm(1, &s[0], &c, &mob[0], &dmob[0]);
const double* mu = props_.viscosity();
mob[0] /= mu[0];
mob[1] /= mu[1];
/// \endcode
/// \page tutorial4
/// \details We
/// recall that we use Fortran ordering for kr derivatives,
/// therefore dmob[i*2 + j] is row j and column i of the
/// matrix.
/// Each row corresponds to a kr function, so which mu to
/// divide by also depends on the row, j.
/// \code
dmob[0*2 + 0] /= mu[0];
dmob[0*2 + 1] /= mu[1];
dmob[1*2 + 0] /= mu[0];
dmob[1*2 + 1] /= mu[1];
}
template <class Sat,
class Pcap,
class DPcap>
void TwophaseFluid::pc(int /*c */, const Sat& /* s*/, Pcap& pcap, DPcap& dpcap) const
{
pcap = 0.;
dpcap = 0.;
}
double TwophaseFluid::s_min(int c) const
{
return smin_[2*c + 0];
}
double TwophaseFluid::s_max(int c) const
{
return smax_[2*c + 0];
}
/// \endcode
/// \page tutorial4 /// \page tutorial4