Loading...
Searching...
No Matches
setupBoundaryConditions.hpp
//===========================================================================
//
// File: setupBoundaryConditions.hpp
//
// Created: Fri Aug 21 09:07:09 2009
//
// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
// Bård Skaflestad <bard.skaflestad@sintef.no>
//
// $Date$
//
// $Revision$
//
//===========================================================================
/*
Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
Copyright 2009, 2010 Statoil ASA.
This file is part of The Open Reservoir Simulator Project (OpenRS).
OpenRS 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.
OpenRS 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 OpenRS. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPENRS_SETUPBOUNDARYCONDITIONS_HEADER
#define OPENRS_SETUPBOUNDARYCONDITIONS_HEADER
#include <opm/common/utility/parameters/ParameterGroup.hpp>
#include <opm/input/eclipse/Units/Units.hpp>
#include <opm/porsol/common/BoundaryConditions.hpp>
#include <opm/porsol/common/PeriodicHelpers.hpp>
namespace Opm
{
template <class GridInterface, class BCs>
inline void setupBoundaryConditions(const Opm::ParameterGroup& param,
const GridInterface& g,
BCs& bcs)
{
if (param.getDefault("upscaling", false)) {
int bct = param.get<int>("boundary_condition_type");
int pddir = param.getDefault("pressure_drop_direction", 0);
double pdrop = param.getDefault("boundary_pressuredrop", 1.0e5);
double bdy_sat = param.getDefault("boundary_saturation", 1.0);
bool twodim_hack = param.getDefault("2d_hack", false);
setupUpscalingConditions(g, bct, pddir, pdrop, bdy_sat, twodim_hack, bcs);
return;
}
if (param.getDefault("region_based_bcs", false)) {
setupRegionBasedConditions(param, g, bcs);
return;
}
// Make flow equation boundary conditions.
// Default is pressure 1.0e5 on the left, 0.0 on the right.
// Recall that the boundary ids range from 1 to 6 for the cartesian edges,
// and that boundary id 0 means interiour face/intersection.
std::string flow_bc_type = param.getDefault<std::string>("flow_bc_type", "dirichlet");
FlowBC::BCType bct = FlowBC::Dirichlet;
double leftval = 1.0*Opm::unit::barsa;
double rightval = 0.0;
if (flow_bc_type == "neumann") {
bct = FlowBC::Neumann;
leftval = param.get<double>("left_flux");
rightval = param.getDefault<double>("right_flux", -leftval);
} else if (flow_bc_type == "dirichlet") {
leftval = param.getDefault<double>("left_pressure", leftval);
rightval = param.getDefault<double>("right_pressure", rightval);
} else if (flow_bc_type == "periodic") {
OPM_THROW(std::runtime_error, "Periodic conditions not here yet.");
} else {
OPM_THROW(std::runtime_error, "Unknown flow boundary condition type " + flow_bc_type);
}
bcs.resize(7);
bcs.flowCond(1) = FlowBC(bct, leftval);
bcs.flowCond(2) = FlowBC(bct, rightval);
// Default transport boundary conditions are used.
}
template <class GridInterface, class BCs>
inline void setupUpscalingConditions(const GridInterface& g,
int bct,
int pddir,
double pdrop,
double bdy_sat,
bool twodim_hack,
BCs& bcs)
{
// Caution: This enum is copied from Upscaler.hpp.
enum BoundaryConditionType { Fixed = 0, Linear = 1, Periodic = 2, PeriodicSingleDirection = 3, Noflow = 4 };
if (bct < 0 || bct > 2) {
OPM_THROW(std::runtime_error,
"Illegal boundary condition type (0-2 are legal): " + std::to_string(bct)); // Later on, we may allow 3 and 4.
}
BoundaryConditionType bctype = static_cast<BoundaryConditionType>(bct);
assert(pddir >=0 && pddir <= 2);
// Flow conditions.
switch (bctype) {
case Fixed:
{
// assert(!g.uniqueBoundaryIds());
bcs.clear();
bcs.resize(7);
bcs.flowCond(2*pddir + 1) = FlowBC(FlowBC::Dirichlet, pdrop);
bcs.flowCond(2*pddir + 2) = FlowBC(FlowBC::Dirichlet, 0.0);
bcs.satCond(2*pddir + 1) = SatBC(SatBC::Dirichlet, bdy_sat); // The only possible inflow location.
for (int i = 0; i < 7; ++i) {
bcs.setCanonicalBoundaryId(i, i);
}
break;
}
case Linear:
{
// assert(g.uniqueBoundaryIds());
createLinear(bcs, g, pdrop, pddir, bdy_sat, twodim_hack);
break;
}
case Periodic:
{
// assert(g.uniqueBoundaryIds());
FlowBC fb(FlowBC::Periodic, 0.0);
std::array<FlowBC, 6> fcond = {{ fb, fb, fb, fb, fb, fb }};
fcond[2*pddir] = FlowBC(FlowBC::Periodic, pdrop);
fcond[2*pddir + 1] = FlowBC(FlowBC::Periodic, -pdrop);
SatBC sb(SatBC::Periodic, 0.0);
std::array<SatBC, 6> scond = {{ sb, sb, sb, sb, sb, sb }};
if (twodim_hack) {
// fcond[2] = FlowBC(FlowBC::Neumann, 0.0);
// fcond[3] = FlowBC(FlowBC::Neumann, 0.0);
fcond[4] = FlowBC(FlowBC::Neumann, 0.0);
fcond[5] = FlowBC(FlowBC::Neumann, 0.0);
// scond[2] = SatBC(SatBC::Dirichlet, 1.0);
// scond[3] = SatBC(SatBC::Dirichlet, 1.0);
scond[4] = SatBC(SatBC::Dirichlet, 1.0);
scond[5] = SatBC(SatBC::Dirichlet, 1.0);
}
createPeriodic(bcs, g, fcond, scond);
break;
}
default:
OPM_THROW(std::runtime_error, "Error in switch statement, should never be here.");
}
// Default transport boundary conditions are used.
}
namespace
{
template <class Vector>
bool isInside(const Vector& low, const Vector& high, const Vector& pt)
{
return low[0] < pt[0]
&& low[1] < pt[1]
&& low[2] < pt[2]
&& high[0] > pt[0]
&& high[1] > pt[1]
&& high[2] > pt[2];
}
} // anon namespace
template <class GridInterface, class BCs>
inline void setupRegionBasedConditions(const Opm::ParameterGroup& param,
const GridInterface& g,
BCs& bcs)
{
// Extract region and pressure value for Dirichlet bcs.
typedef typename GridInterface::Vector Vector;
Vector low;
low[0] = param.getDefault("dir_block_low_x", 0.0);
low[1] = param.getDefault("dir_block_low_y", 0.0);
low[2] = param.getDefault("dir_block_low_z", 0.0);
Vector high;
high[0] = param.getDefault("dir_block_high_x", 1.0);
high[1] = param.getDefault("dir_block_high_y", 1.0);
high[2] = param.getDefault("dir_block_high_z", 1.0);
double dir_block_pressure = param.get<double>("dir_block_pressure");
// Set flow conditions for that region.
// For this to work correctly, unique boundary ids should be used,
// otherwise conditions may spread outside the given region, to all
// faces with the same bid as faces inside the region.
typedef typename GridInterface::CellIterator CI;
typedef typename CI::FaceIterator FI;
int max_bid = 0;
std::vector<int> dir_bids;
for (CI c = g.cellbegin(); c != g.cellend(); ++c) {
for (FI f = c->facebegin(); f != c->faceend(); ++f) {
int bid = f->boundaryId();
max_bid = std::max(bid, max_bid);
if (bid != 0 && isInside(low, high, f->centroid())) {
dir_bids.push_back(bid);
}
}
}
bcs.resize(max_bid + 1);
for (std::vector<int>::const_iterator it = dir_bids.begin(); it != dir_bids.end(); ++it) {
bcs.flowCond(*it) = FlowBC(FlowBC::Dirichlet, dir_block_pressure);
}
// Transport BCs are defaulted.
}
} // namespace Opm
#endif // OPENRS_SETUPBOUNDARYCONDITIONS_HEADER
Opm::BCBase< double >::BCType
BCType
Definition BoundaryConditions.hpp:58
Opm::FlowBC
A class for representing a flow boundary condition.
Definition BoundaryConditions.hpp:122
Opm::SatBC
A class for representing a saturation boundary condition.
Definition BoundaryConditions.hpp:176
Opm::GridInterface
Opm::Elasticity::Vector
Dune::BlockVector< Dune::FieldVector< double, 1 > > Vector
A vector holding our RHS.
Definition matrixops.hpp:33
Opm
Inverting small matrices.
Definition ImplicitAssembly.hpp:43
Opm::createLinear
void createLinear(BCs &fbcs, const GridInterface &g, const double pdrop, const int pddir, const double bdy_sat, const bool twodim_hack=false, const double spatial_tolerance=1e-6)
Makes a boundary condition object representing linear boundary conditions in any cartesian direction.
Definition PeriodicHelpers.hpp:220
Opm::setupBoundaryConditions
void setupBoundaryConditions(const Opm::ParameterGroup &param, const GridInterface &g, BCs &bcs)
Setup boundary conditions for a simulation.
Definition setupBoundaryConditions.hpp:51
Opm::setupUpscalingConditions
void setupUpscalingConditions(const GridInterface &g, int bct, int pddir, double pdrop, double bdy_sat, bool twodim_hack, BCs &bcs)
Definition setupBoundaryConditions.hpp:99
Opm::createPeriodic
void createPeriodic(BCs &fbcs, const GridInterface &g, const std::array< FlowBC, 2 *GridInterface::Dimension > &fconditions, const std::array< SatBC, 2 *GridInterface::Dimension > &sconditions, double spatial_tolerance=1e-6)
Makes a boundary condition object representing periodic boundary conditions in any cartesian directio...
Definition PeriodicHelpers.hpp:105