setupBoundaryConditions.hpp
Go to the documentation of this file.
1//===========================================================================
2//
3// File: setupBoundaryConditions.hpp
4//
5// Created: Fri Aug 21 09:07:09 2009
6//
7// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
8// Bård Skaflestad <bard.skaflestad@sintef.no>
9//
10// $Date$
11//
12// $Revision$
13//
14//===========================================================================
15
16/*
17 Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
18 Copyright 2009, 2010 Statoil ASA.
19
20 This file is part of The Open Reservoir Simulator Project (OpenRS).
21
22 OpenRS is free software: you can redistribute it and/or modify
23 it under the terms of the GNU General Public License as published by
24 the Free Software Foundation, either version 3 of the License, or
25 (at your option) any later version.
26
27 OpenRS is distributed in the hope that it will be useful,
28 but WITHOUT ANY WARRANTY; without even the implied warranty of
29 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 GNU General Public License for more details.
31
32 You should have received a copy of the GNU General Public License
33 along with OpenRS. If not, see <http://www.gnu.org/licenses/>.
34*/
35
36#ifndef OPENRS_SETUPBOUNDARYCONDITIONS_HEADER
37#define OPENRS_SETUPBOUNDARYCONDITIONS_HEADER
38
39#include <opm/core/utility/parameters/ParameterGroup.hpp>
40#include <opm/core/utility/Units.hpp>
41#include <opm/porsol/common/BoundaryConditions.hpp>
42#include <opm/porsol/common/PeriodicHelpers.hpp>
43
44namespace Opm
45{
46
50 template <class GridInterface, class BCs>
51 inline void setupBoundaryConditions(const Opm::parameter::ParameterGroup& param,
52 const GridInterface& g,
53 BCs& bcs)
54 {
55 if (param.getDefault("upscaling", false)) {
56 int bct = param.get<int>("boundary_condition_type");
57 int pddir = param.getDefault("pressure_drop_direction", 0);
58 double pdrop = param.getDefault("boundary_pressuredrop", 1.0e5);
59 double bdy_sat = param.getDefault("boundary_saturation", 1.0);
60 bool twodim_hack = param.getDefault("2d_hack", false);
61 setupUpscalingConditions(g, bct, pddir, pdrop, bdy_sat, twodim_hack, bcs);
62 return;
63 }
64 if (param.getDefault("region_based_bcs", false)) {
65 setupRegionBasedConditions(param, g, bcs);
66 return;
67 }
68 // Make flow equation boundary conditions.
69 // Default is pressure 1.0e5 on the left, 0.0 on the right.
70 // Recall that the boundary ids range from 1 to 6 for the cartesian edges,
71 // and that boundary id 0 means interiour face/intersection.
72 std::string flow_bc_type = param.getDefault<std::string>("flow_bc_type", "dirichlet");
73 FlowBC::BCType bct = FlowBC::Dirichlet;
74 double leftval = 1.0*Opm::unit::barsa;
75 double rightval = 0.0;
76 if (flow_bc_type == "neumann") {
77 bct = FlowBC::Neumann;
78 leftval = param.get<double>("left_flux");
79 rightval = param.getDefault<double>("right_flux", -leftval);
80 } else if (flow_bc_type == "dirichlet") {
81 leftval = param.getDefault<double>("left_pressure", leftval);
82 rightval = param.getDefault<double>("right_pressure", rightval);
83 } else if (flow_bc_type == "periodic") {
84 OPM_THROW(std::runtime_error, "Periodic conditions not here yet.");
85 } else {
86 OPM_THROW(std::runtime_error, "Unknown flow boundary condition type " << flow_bc_type);
87 }
88 bcs.resize(7);
89 bcs.flowCond(1) = FlowBC(bct, leftval);
90 bcs.flowCond(2) = FlowBC(bct, rightval);
91
92 // Default transport boundary conditions are used.
93 }
94
98 template <class GridInterface, class BCs>
99 inline void setupUpscalingConditions(const GridInterface& g,
100 int bct,
101 int pddir,
102 double pdrop,
103 double bdy_sat,
104 bool twodim_hack,
105 BCs& bcs)
106 {
107 // Caution: This enum is copied from Upscaler.hpp.
108 enum BoundaryConditionType { Fixed = 0, Linear = 1, Periodic = 2, PeriodicSingleDirection = 3, Noflow = 4 };
109 if (bct < 0 || bct > 2) {
110 OPM_THROW(std::runtime_error, "Illegal boundary condition type (0-2 are legal): " << bct); // Later on, we may allow 3 and 4.
111 }
112 BoundaryConditionType bctype = static_cast<BoundaryConditionType>(bct);
113 assert(pddir >=0 && pddir <= 2);
114
115 // Flow conditions.
116 switch (bctype) {
117 case Fixed:
118 {
119 // assert(!g.uniqueBoundaryIds());
120 bcs.clear();
121 bcs.resize(7);
122 bcs.flowCond(2*pddir + 1) = FlowBC(FlowBC::Dirichlet, pdrop);
123 bcs.flowCond(2*pddir + 2) = FlowBC(FlowBC::Dirichlet, 0.0);
124 bcs.satCond(2*pddir + 1) = SatBC(SatBC::Dirichlet, bdy_sat); // The only possible inflow location.
125 for (int i = 0; i < 7; ++i) {
126 bcs.setCanonicalBoundaryId(i, i);
127 }
128 break;
129 }
130 case Linear:
131 {
132 // assert(g.uniqueBoundaryIds());
133 createLinear(bcs, g, pdrop, pddir, bdy_sat, twodim_hack);
134 break;
135 }
136 case Periodic:
137 {
138 // assert(g.uniqueBoundaryIds());
139 FlowBC fb(FlowBC::Periodic, 0.0);
140 std::array<FlowBC, 6> fcond = {{ fb, fb, fb, fb, fb, fb }};
141 fcond[2*pddir] = FlowBC(FlowBC::Periodic, pdrop);
142 fcond[2*pddir + 1] = FlowBC(FlowBC::Periodic, -pdrop);
143 SatBC sb(SatBC::Periodic, 0.0);
144 std::array<SatBC, 6> scond = {{ sb, sb, sb, sb, sb, sb }};
145 if (twodim_hack) {
146// fcond[2] = FlowBC(FlowBC::Neumann, 0.0);
147// fcond[3] = FlowBC(FlowBC::Neumann, 0.0);
148 fcond[4] = FlowBC(FlowBC::Neumann, 0.0);
149 fcond[5] = FlowBC(FlowBC::Neumann, 0.0);
150// scond[2] = SatBC(SatBC::Dirichlet, 1.0);
151// scond[3] = SatBC(SatBC::Dirichlet, 1.0);
152 scond[4] = SatBC(SatBC::Dirichlet, 1.0);
153 scond[5] = SatBC(SatBC::Dirichlet, 1.0);
154 }
155 createPeriodic(bcs, g, fcond, scond);
156 break;
157 }
158 default:
159 OPM_THROW(std::runtime_error, "Error in switch statement, should never be here.");
160 }
161
162 // Default transport boundary conditions are used.
163 }
164
165
166
167 namespace
168 {
169 template <class Vector>
170 bool isInside(const Vector& low, const Vector& high, const Vector& pt)
171 {
172 return low[0] < pt[0]
173 && low[1] < pt[1]
174 && low[2] < pt[2]
175 && high[0] > pt[0]
176 && high[1] > pt[1]
177 && high[2] > pt[2];
178 }
179 } // anon namespace
180
181
182 template <class GridInterface, class BCs>
183 inline void setupRegionBasedConditions(const Opm::parameter::ParameterGroup& param,
184 const GridInterface& g,
185 BCs& bcs)
186 {
187 // Extract region and pressure value for Dirichlet bcs.
188 typedef typename GridInterface::Vector Vector;
189 Vector low;
190 low[0] = param.getDefault("dir_block_low_x", 0.0);
191 low[1] = param.getDefault("dir_block_low_y", 0.0);
192 low[2] = param.getDefault("dir_block_low_z", 0.0);
193 Vector high;
194 high[0] = param.getDefault("dir_block_high_x", 1.0);
195 high[1] = param.getDefault("dir_block_high_y", 1.0);
196 high[2] = param.getDefault("dir_block_high_z", 1.0);
197 double dir_block_pressure = param.get<double>("dir_block_pressure");
198
199 // Set flow conditions for that region.
200 // For this to work correctly, unique boundary ids should be used,
201 // otherwise conditions may spread outside the given region, to all
202 // faces with the same bid as faces inside the region.
203 typedef typename GridInterface::CellIterator CI;
204 typedef typename CI::FaceIterator FI;
205 int max_bid = 0;
206 std::vector<int> dir_bids;
207 for (CI c = g.cellbegin(); c != g.cellend(); ++c) {
208 for (FI f = c->facebegin(); f != c->faceend(); ++f) {
209 int bid = f->boundaryId();
210 max_bid = std::max(bid, max_bid);
211 if (bid != 0 && isInside(low, high, f->centroid())) {
212 dir_bids.push_back(bid);
213 }
214 }
215 }
216 bcs.resize(max_bid + 1);
217 for (std::vector<int>::const_iterator it = dir_bids.begin(); it != dir_bids.end(); ++it) {
218 bcs.flowCond(*it) = FlowBC(FlowBC::Dirichlet, dir_block_pressure);
219 }
220
221 // Transport BCs are defaulted.
222 }
223
224
225
226} // namespace Opm
227
228
229#endif // OPENRS_SETUPBOUNDARYCONDITIONS_HEADER
Opm::BCBase< double >::BCType
BCType
Enum for the allowed boundary condition types. So far, we support Dirichlet, Neumann and Periodic con...
Definition: BoundaryConditions.hpp:58
Opm::BCBase< double >::Periodic
@ Periodic
Definition: BoundaryConditions.hpp:58
Opm::BCBase< double >::Neumann
@ Neumann
Definition: BoundaryConditions.hpp:58
Opm::BCBase< double >::Dirichlet
@ Dirichlet
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
Definition: BlackoilFluid.hpp:32
Opm::setupBoundaryConditions
void setupBoundaryConditions(const Opm::parameter::ParameterGroup &param, const GridInterface &g, BCs &bcs)
Setup boundary conditions for a simulation. It is assumed that the boundary ids are 1-6,...
Definition: setupBoundaryConditions.hpp:51
Opm::setupRegionBasedConditions
void setupRegionBasedConditions(const Opm::parameter::ParameterGroup &param, const GridInterface &g, BCs &bcs)
Definition: setupBoundaryConditions.hpp:183
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:224
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