diff --git a/generate_doc_figures.py b/generate_doc_figures.py index 8e53d167..ae42709d 100755 --- a/generate_doc_figures.py +++ b/generate_doc_figures.py @@ -90,6 +90,34 @@ for case in cases: WriteImage(join(figure_path, "tutorial3-"+case+".png")) 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 for f in collected_garbage_file: remove(f) diff --git a/opm/core/grid.h b/opm/core/grid.h index b6dfb3f0..841ea8e2 100644 --- a/opm/core/grid.h +++ b/opm/core/grid.h @@ -250,18 +250,20 @@ struct UnstructuredGrid * create_grid_empty(void); /** - Allocate and initialise an UnstructuredGrid where pointers are set to - location with correct size. + Allocate and initialise an UnstructuredGrid where pointers are set + to location with correct size. \param[in] ndims Number of physical dimensions. \param[in] ncells Number of cells. \param[in] nfaces Number of faces. \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] nnodes Number of Nodes - \return Fully formed UnstructuredGrid with all fields allocated, but not filled with - usefull numbers. - NULL in case of allocation failure. + \param[in] ncellfaces Size of mapping from cells to faces. + (i.e., the number of `half-faces') + \param[in] nnodes Number of Nodes. + + \return Fully formed UnstructuredGrid with all fields except + global_cell allocated, but not filled with meaningful + values. NULL in case of allocation failure. */ struct UnstructuredGrid * allocate_grid(size_t ndims , diff --git a/opm/core/grid/cart_grid.c b/opm/core/grid/cart_grid.c index b318de74..2847e1fd 100644 --- a/opm/core/grid/cart_grid.c +++ b/opm/core/grid/cart_grid.c @@ -170,68 +170,13 @@ allocate_cart_grid(size_t ndims , size_t nfaces, size_t nnodes) { - size_t nel; - struct UnstructuredGrid *G; + size_t nfacenodes, ncellfaces; - G = create_grid_empty(); + nfacenodes = nfaces * (2 * (ndims - 1)); + ncellfaces = ncells * (2 * ndims); - if (G != NULL) { - /* Node fields ---------------------------------------- */ - 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; + return allocate_grid(ndims, ncells, nfaces, + nfacenodes, ncellfaces, nnodes); } diff --git a/opm/core/linalg/LinearSolverFactory.cpp b/opm/core/linalg/LinearSolverFactory.cpp index 6943b88c..9b6d8909 100644 --- a/opm/core/linalg/LinearSolverFactory.cpp +++ b/opm/core/linalg/LinearSolverFactory.cpp @@ -55,20 +55,29 @@ namespace Opm LinearSolverFactory::LinearSolverFactory(const parameter::ParameterGroup& param) { - const std::string ls = param.getDefault("linsolver", "umfpack"); + const std::string ls = + param.getDefault("linsolver", "umfpack"); + if (ls == "umfpack") { #if HAVE_SUITESPARSE_UMFPACK_H solver_.reset(new LinearSolverUmfpack); -#else - THROW("Linear solver " << ls <<" is not available."); #endif - } else if (ls == "istl") { + } + + else if (ls == "istl") { #if HAVE_DUNE_ISTL solver_.reset(new LinearSolverIstl(param)); -#else - THROW("Linear solver " << ls <<" is not available."); #endif } + + else { + THROW("Linear solver " << ls << " is unknown."); + } + + if (! solver_) { + THROW("Linear solver " << ls << " is not enabled in " + "this configuration."); + } } diff --git a/opm/core/wells/WellsGroup.cpp b/opm/core/wells/WellsGroup.cpp index 3514bb5d..3da16f77 100644 --- a/opm/core/wells/WellsGroup.cpp +++ b/opm/core/wells/WellsGroup.cpp @@ -941,7 +941,6 @@ namespace Opm double WellNode::productionGuideRate(bool only_group) { if (!only_group || prodSpec().control_mode_ == ProductionSpecification::GRUP) { - std::cout << prodSpec().guide_rate_ << std::endl; return prodSpec().guide_rate_; } return 0.0; diff --git a/tutorials/tutorial4.cpp b/tutorials/tutorial4.cpp index 7bdbd1b9..9358da13 100644 --- a/tutorials/tutorial4.cpp +++ b/tutorials/tutorial4.cpp @@ -35,14 +35,7 @@ #include #include -#include -#include -#include -#include -#include -#include -#include -#include +#include #include @@ -52,113 +45,6 @@ #include /// \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 - 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 - 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 smin_; - std::vector smax_; -}; -/// \endcode - -/// \page tutorial4 -/// \details We set up the transport model. -/// \code -typedef Opm::SinglePointUpwindTwoPhase 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 > NVecColl; -typedef JacobianSystem< struct CSRMatrix, NVecColl > JacSys; - -template -class MaxNorm { -public: - static double - norm(const Vector& v) { - return AccumulationNorm ::norm(v); - } -}; -/// \endcode - -/// \page tutorial4 -/// \details -/// We set up the solver. -/// \code -typedef Opm::ImplicitTransport TransportSolver; -/// \endcode - /// \page tutorial4 /// \details @@ -179,8 +65,9 @@ int main () double dy = 10.; double dz = 10.; using namespace Opm; - GridManager grid(nx, ny, nz, dx, dy, dz); - int num_cells = grid.c_grid()->number_of_cells; + GridManager grid_manager(nx, ny, nz, dx, dy, dz); + const UnstructuredGrid& grid = *grid_manager.c_grid(); + int num_cells = grid.number_of_cells; /// \endcode /// \page tutorial4 @@ -214,8 +101,7 @@ int main () /// description and set this function to be linear. For more realistic fluid, the /// saturation function is given by the data. /// \code - SaturationPropsBasic::RelPermFunc rel_perm_func; - rel_perm_func = SaturationPropsBasic::Linear; + SaturationPropsBasic::RelPermFunc rel_perm_func = SaturationPropsBasic::Linear; /// \endcode /// \page tutorial4 @@ -224,7 +110,6 @@ int main () /// \code IncompPropertiesBasic props(num_phases, rel_perm_func, rho, mu, porosity, k, dim, num_cells); - TwophaseFluid fluid(props); /// \endcode @@ -236,15 +121,6 @@ int main () /// \endcode - - - /// \page tutorial4 - /// \details We set up the source term - /// \code - std::vector src(num_cells, 0.0); - src[0] = 1.; - src[num_cells-1] = -1.; - /// \endcode /// \page tutorial4 /// \details We set up necessary information for the wells @@ -259,34 +135,27 @@ int main () /// \endcode /// \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 - TransportSource* tsrc = create_transport_source(2, 2); - double ssrc[] = { 1.0, 0.0 }; - double ssink[] = { 0.0, 1.0 }; - 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); - } - } + std::vector src(num_cells, 0.0); + src[0] = 1.; + src[num_cells-1] = -1.; /// \endcode /// \page tutorial4 /// \details We compute the pore volume /// \code std::vector porevol; - Opm::computePorevolume(*grid.c_grid(), props.porosity(), porevol); + Opm::computePorevolume(grid, props.porosity(), porevol); /// \endcode /// \page tutorial4 - /// \details We set up the transport solver. + /// \details Set up the transport solver. This is a reordering implicit Euler transport solver. /// \code - const bool guess_old_solution = true; - TransportModel model (fluid, *grid.c_grid(), porevol, grav, guess_old_solution); - TransportSolver tsolver(model); + const double tolerance = 1e-9; + const int max_iterations = 30; + Opm::TransportModelTwophase transport_solver(grid, props, tolerance, max_iterations); /// \endcode /// \page tutorial4 @@ -296,14 +165,6 @@ int main () int num_time_steps = 20; /// \endcode - /// \page tutorial4 - /// \details Control paramaters for the implicit solver. - /// \code - ImplicitTransportDetails::NRReport rpt; - ImplicitTransportDetails::NRControl ctrl; - /// \endcode - - /// \page tutorial4 /// \details We define a vector which contains all cell indexes. We use this @@ -321,21 +182,13 @@ int main () FlowBCManager bcs; /// \endcode - /// \page tutorial4 - /// \details - /// Linear solver init. - /// \code - using ImplicitTransportLinAlgSupport::CSRMatrixUmfpackSolver; - CSRMatrixUmfpackSolver linsolve; - /// \endcode - /// \page tutorial4 /// \details /// We set up a two-phase state object, and /// initialise water saturation to minimum everywhere. /// \code TwophaseState state; - state.init(*grid.c_grid(), 2); + state.init(grid, 2); state.setFirstSat(allcells, props, TwophaseState::MinSat); /// \endcode @@ -450,7 +303,7 @@ int main () /// last argument. /// \code LinearSolverUmfpack linsolver; - IncompTpfa psolver(*grid.c_grid(), props.permeability(), grav, linsolver, wells); + IncompTpfa psolver(grid, props.permeability(), grav, linsolver, wells); /// \endcode @@ -469,7 +322,7 @@ int main () /// \page tutorial4 /// \details In order to use the well controls, we need to generate the WDP for each well. /// \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 /// \page tutorial4 @@ -511,7 +364,7 @@ int main () /// \page tutorial4 /// \details Transport solver /// \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 /// \page tutorial4 @@ -523,75 +376,11 @@ int main () Opm::DataMap dm; dm["saturation"] = &state.saturation(); dm["pressure"] = &state.pressure(); - Opm::writeVtkData(*grid.c_grid(), dm, vtkfile); + Opm::writeVtkData(grid, dm, vtkfile); } } /// \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 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 -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 -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