240 lines
7.6 KiB
C++
240 lines
7.6 KiB
C++
/*
|
|
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/>.
|
|
*/
|
|
|
|
|
|
#include "config.h"
|
|
#include <opm/core/pressure/FlowBCManager.hpp>
|
|
#include <opm/common/ErrorMacros.hpp>
|
|
#include <opm/core/grid.h>
|
|
#include <vector>
|
|
|
|
namespace Opm
|
|
{
|
|
|
|
namespace
|
|
{
|
|
std::string sideString(FlowBCManager::Side s);
|
|
void findSideFaces(const UnstructuredGrid& grid,
|
|
const FlowBCManager::Side side,
|
|
std::vector<int>& faces);
|
|
} // anon namespace
|
|
|
|
|
|
/// Default constructor sets up empty boundary conditions.
|
|
/// By convention, this is equivalent to all-noflow conditions.
|
|
FlowBCManager::FlowBCManager()
|
|
: bc_(0)
|
|
{
|
|
bc_ = flow_conditions_construct(0);
|
|
if (!bc_) {
|
|
OPM_THROW(std::runtime_error, "Failed to construct FlowBoundaryConditions struct.");
|
|
}
|
|
}
|
|
|
|
|
|
/// Destructor.
|
|
FlowBCManager::~FlowBCManager()
|
|
{
|
|
flow_conditions_destroy(bc_);
|
|
}
|
|
|
|
|
|
/// Remove all appended BCs.
|
|
/// By convention, BCs are now equivalent to all-noflow conditions.
|
|
void FlowBCManager::clear()
|
|
{
|
|
flow_conditions_clear(bc_);
|
|
}
|
|
|
|
|
|
/// Append a single boundary condition.
|
|
/// If the type is BC_NOFLOW the value argument is not used.
|
|
/// If the type is BC_PRESSURE the value argument is a pressure value.
|
|
/// If the type is BC_FLUX_TOTVOL the value argument is a total flux value (m^3/s).
|
|
/// Note: unset boundary conditions are noflow by convention,
|
|
/// so it is normally not necessary to explicitly append
|
|
/// BC_NOFLOW conditions. However, it may make sense to do so
|
|
/// if the bc will change during a simulation run.
|
|
/// Note: if normal velocity bcs are desired, convert to
|
|
/// fluxes by multiplying with face area.
|
|
void FlowBCManager::append(const FlowBCType type,
|
|
const int face,
|
|
const double value)
|
|
{
|
|
int ok = flow_conditions_append(type, face, value, bc_);
|
|
if (!ok) {
|
|
OPM_THROW(std::runtime_error, "Failed to append boundary condition for face " << face);
|
|
}
|
|
}
|
|
|
|
|
|
/// Add BC_PRESSURE boundary conditions to all faces on a given side.
|
|
/// The grid must have a logical cartesian structure, and grid
|
|
/// faces must be tagged (i.e. grid.cell_facetag must be
|
|
/// non-null). Only the set of faces adjacent to cells with
|
|
/// minimum/maximum I/J/K coordinate (depending on side) are
|
|
/// considered.
|
|
void FlowBCManager::pressureSide(const UnstructuredGrid& grid,
|
|
const Side side,
|
|
const double pressure)
|
|
{
|
|
std::vector<int> faces;
|
|
findSideFaces(grid, side, faces);
|
|
int ok = flow_conditions_append_multi(BC_PRESSURE, faces.size(), &faces[0], pressure, bc_);
|
|
if (!ok) {
|
|
OPM_THROW(std::runtime_error, "Failed to append pressure boundary conditions for side " << sideString(side));
|
|
}
|
|
}
|
|
|
|
|
|
/// Add BC_FLUX_TOTVOL boundary conditions to all faces on a given side.
|
|
/// The grid must have a logical cartesian structure, and grid
|
|
/// faces must be tagged (i.e. grid.cell_facetag must be
|
|
/// non-null). Only the set of faces adjacent to cells with
|
|
/// minimum/maximum I/J/K coordinate (depending on side) are
|
|
/// considered.
|
|
/// The flux specified is taken to be the total flux through
|
|
/// the side, each individual face receiving a part of the
|
|
/// total flux in proportion to its area, so that all faces
|
|
/// will have identical normal velocities.
|
|
void FlowBCManager::fluxSide(const UnstructuredGrid& grid,
|
|
const Side side,
|
|
const double flux)
|
|
{
|
|
// Find side faces.
|
|
std::vector<int> faces;
|
|
findSideFaces(grid, side, faces);
|
|
|
|
// Compute total area of faces.
|
|
double tot_area = 0.0;
|
|
for (int fi = 0; fi < int(faces.size()); ++fi) {
|
|
tot_area += grid.face_areas[faces[fi]];
|
|
}
|
|
|
|
// Append flux conditions for all the faces individually.
|
|
for (int fi = 0; fi < int(faces.size()); ++fi) {
|
|
const double face_flux = flux * grid.face_areas[faces[fi]] / tot_area;
|
|
int ok = flow_conditions_append(BC_FLUX_TOTVOL, faces[fi], face_flux, bc_);
|
|
if (!ok) {
|
|
OPM_THROW(std::runtime_error, "Failed to append flux boundary conditions for face " << faces[fi] << " on side " << sideString(side));
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
|
|
/// Access the managed boundary conditions.
|
|
/// The method is named similarly to c_str() in std::string,
|
|
/// to make it clear that we are returning a C-compatible struct.
|
|
const FlowBoundaryConditions* FlowBCManager::c_bcs() const
|
|
{
|
|
return bc_;
|
|
}
|
|
|
|
|
|
|
|
|
|
// ------ Utility functions ------
|
|
|
|
|
|
namespace
|
|
{
|
|
std::string sideString(FlowBCManager::Side s)
|
|
{
|
|
switch (s) {
|
|
case FlowBCManager::Xmin: return "Xmin";
|
|
case FlowBCManager::Xmax: return "Xmax";
|
|
case FlowBCManager::Ymin: return "Ymin";
|
|
case FlowBCManager::Ymax: return "Ymax";
|
|
case FlowBCManager::Zmin: return "Zmin";
|
|
case FlowBCManager::Zmax: return "Zmax";
|
|
default: OPM_THROW(std::runtime_error, "Unknown side tag " << s);
|
|
}
|
|
}
|
|
|
|
|
|
|
|
void cartCoord(const int ndims,
|
|
const int log_cart_coord,
|
|
const int* dims,
|
|
int* ijk)
|
|
{
|
|
int ix = log_cart_coord;
|
|
|
|
for (int dim = 0; dim < ndims; ++dim) {
|
|
ijk[dim] = ix % dims[dim];
|
|
ix /= dims[dim];
|
|
}
|
|
|
|
// Make sure that lexicographic index is consistent with
|
|
// grid dimensions.
|
|
assert(ix == 0);
|
|
}
|
|
|
|
|
|
|
|
/// The grid must have a logical cartesian structure, and grid
|
|
/// faces must be tagged (i.e. grid.cell_facetag must be
|
|
/// non-null). Only the set of faces adjacent to cells with
|
|
/// minimum/maximum I/J/K coordinate (depending on side) are
|
|
/// considered.
|
|
void findSideFaces(const UnstructuredGrid& grid,
|
|
const FlowBCManager::Side side,
|
|
std::vector<int>& faces)
|
|
{
|
|
if (grid.cell_facetag == 0) {
|
|
OPM_THROW(std::runtime_error, "Faces not tagged - cannot extract " << sideString(side) << " faces.");
|
|
}
|
|
|
|
// make sure that grid has three dimensions or less.
|
|
assert(grid.dimensions <= 3);
|
|
|
|
// Make sure boundary condition side is consistent with
|
|
// number of physical grid dimensions.
|
|
assert(side < 2 * grid.dimensions);
|
|
|
|
// Get all boundary faces with the correct tag and with
|
|
// min/max i/j/k (depending on side).
|
|
const int correct_ijk = (side % 2) ? grid.cartdims[side/2] - 1 : 0;
|
|
for (int c = 0; c < grid.number_of_cells; ++c) {
|
|
int ijk[3] = { -1, -1, -1 };
|
|
int gc = (grid.global_cell != 0) ? grid.global_cell[c] : c;
|
|
cartCoord(grid.dimensions, gc, grid.cartdims, ijk);
|
|
if (ijk[side/2] != correct_ijk) {
|
|
continue;
|
|
}
|
|
for (int hf = grid.cell_facepos[c]; hf < grid.cell_facepos[c + 1]; ++hf) {
|
|
if (grid.cell_facetag[hf] == side) {
|
|
// Tag is correct.
|
|
const int f = grid.cell_faces[hf];
|
|
if (grid.face_cells[2*f] == -1 || grid.face_cells[2*f + 1] == -1) {
|
|
// Face is on boundary.
|
|
faces.push_back(f);
|
|
} else {
|
|
OPM_THROW(std::runtime_error, "Face not on boundary, even with correct tag and boundary cell. This should not occur.");
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
} // anon namespace
|
|
|
|
} // namespace Opm
|