36#ifndef OPENRS_PERIODICHELPERS_HEADER
37#define OPENRS_PERIODICHELPERS_HEADER
41#include <dune/common/fvector.hh>
55 const std::vector<BoundaryFaceInfo>& bfinfo,
56 const std::array<FlowBC, 6>& fconditions,
57 const std::array<double, 6>& side_areas)
59 int num_bdy = bfinfo.size();
60 for (
int i = 0; i < num_bdy; ++i) {
61 FlowBC face_bc = fconditions[bfinfo[i].canon_pos];
65 bcs.flowCond(bfinfo[i].bid) = face_bc;
75 const std::vector<BoundaryFaceInfo>& bfinfo,
76 const std::array<SatBC, 6>& sconditions)
78 int num_bdy = bfinfo.size();
79 for (
int i = 0; i < num_bdy; ++i) {
80 bcs.satCond(bfinfo[i].bid) = sconditions[bfinfo[i].canon_pos];
88 std::array<bool, 6> retval = {{ bcs[0].isPeriodic(),
93 bcs[5].isPeriodic() }};
104 template <
class BCs,
class Gr
idInterface>
106 const GridInterface& g,
107 const std::array<FlowBC, 2*GridInterface::Dimension>& fconditions,
108 const std::array<SatBC, 2*GridInterface::Dimension>& sconditions,
109 double spatial_tolerance = 1e-6)
111 static_assert(BCs::HasFlowConds,
"");
112 static_assert(BCs::HasSatConds,
"");
114 for (
int i = 0; i < GridInterface::Dimension; ++i) {
115 if (fconditions[2*i].isPeriodic()) {
116 assert(fconditions[2*i + 1].isPeriodic());
117 assert(fconditions[2*i].pressureDifference() == -fconditions[2*i + 1].pressureDifference());
118 assert(sconditions[2*i].isPeriodic());
119 assert(sconditions[2*i + 1].isPeriodic());
120 assert(sconditions[2*i].saturationDifference() == 0.0);
121 assert(sconditions[2*i + 1].saturationDifference() == 0.0);
124 std::vector<BoundaryFaceInfo> bfinfo;
125 std::array<double, 6> side_areas;
134 template <
class BCs,
class Gr
idInterface>
136 const GridInterface& g,
137 const std::array<FlowBC, 2*GridInterface::Dimension>& fconditions,
138 double spatial_tolerance = 1e-6)
140 static_assert(BCs::HasFlowConds,
"");
141 static_assert(!BCs::HasSatConds,
"");
143 for (
int i = 0; i < GridInterface::Dimension; ++i) {
144 if (fconditions[2*i].isPeriodic()) {
145 assert(fconditions[2*i + 1].isPeriodic());
146 assert(fconditions[2*i].pressureDifference() == -fconditions[2*i + 1].pressureDifference());
149 std::vector<BoundaryFaceInfo> bfinfo;
150 std::array<double, 6> side_areas;
158 template <
class BCs,
class Gr
idInterface>
160 const GridInterface& g,
161 const std::array<SatBC, 2*GridInterface::Dimension>& sconditions,
162 double spatial_tolerance = 1e-6)
164 static_assert(!BCs::HasFlowConds,
"");
165 static_assert(BCs::HasSatConds,
"");
167 for (
int i = 0; i < GridInterface::Dimension; ++i) {
168 if (sconditions[2*i].isPeriodic()) {
169 assert(sconditions[2*i + 1].isPeriodic());
170 assert(sconditions[2*i].saturationDifference() == -sconditions[2*i + 1].saturationDifference());
173 std::vector<BoundaryFaceInfo> bfinfo;
174 std::array<double, 6> side_areas;
182 template <
class BCs,
class Gr
idInterface>
184 std::vector<BoundaryFaceInfo>& bfinfo,
185 std::array<double, 6>& side_areas,
186 const GridInterface& g,
187 const std::array<bool, 2*GridInterface::Dimension>& is_periodic,
188 double spatial_tolerance = 1e-6)
190#if DUNE_VERSION_NEWER(DUNE_GRID, 2, 3)
191 findPeriodicPartners(bfinfo, side_areas, g.grid().leafGridView(), is_periodic, spatial_tolerance);
195 int num_bdy = bfinfo.size();
197 int max_bid = num_bdy;
201 fbcs.resize(max_bid + 1);
204 for (
int i = 0; i < num_bdy; ++i) {
205 int bid1 = bfinfo[i].bid;
206 int bid2 = bfinfo[i].partner_bid;
209 fbcs.setPeriodicPartners(bid1, bid2);
211 fbcs.setCanonicalBoundaryId(bid1, bfinfo[i].canon_pos + 1);
223 template <
class BCs,
class Gr
idInterface>
225 const GridInterface& g,
228 const double bdy_sat,
229 const bool twodim_hack =
false,
230 const double spatial_tolerance = 1e-6)
235 typedef typename GridInterface::CellIterator CI;
236 typedef typename CI::FaceIterator FI;
237 typedef typename GridInterface::Vector Vector;
238 std::vector<FI> bface_iters;
242 for (CI c = g.cellbegin(); c != g.cellend(); ++c) {
243 for (FI f = c->facebegin(); f != c->faceend(); ++f) {
244 if (f->boundaryId()) {
245 bface_iters.push_back(f);
246 Vector fcent = f->centroid();
247 for (
int dd = 0; dd < GridInterface::Dimension; ++dd) {
248 low[dd] = std::min(low[dd], fcent[dd]);
249 hi[dd] = std::max(hi[dd], fcent[dd]);
251 max_bid = std::max(max_bid, f->boundaryId());
255 int num_bdy = bface_iters.size();
256 if (max_bid != num_bdy) {
257 OPM_THROW(std::runtime_error,
"createLinear() assumes that every boundary face has a unique boundary id. That seems to be violated.");
259 fbcs.resize(max_bid + 1);
262 double cmin = low[pddir];
263 double cmax = hi[pddir];
264 double cdelta = cmax - cmin;
266 for (
int i = 0; i < num_bdy; ++i) {
267 Vector fcent = bface_iters[i]->centroid();
269 for (
int dd = 0; dd < GridInterface::Dimension; ++dd) {
270 double coord = fcent[dd];
271 if (fabs(coord - low[dd]) <= spatial_tolerance) {
274 }
else if (fabs(coord - hi[dd]) <= spatial_tolerance) {
275 canon_pos = 2*dd + 1;
279 if (canon_pos == -1) {
280 std::cerr <<
"Centroid: " << fcent <<
"\n";
281 std::cerr <<
"Bounding box min: " << low <<
"\n";
282 std::cerr <<
"Bounding box max: " << hi <<
"\n";
283 OPM_THROW(std::runtime_error,
"Boundary face centroid not on bounding box. Maybe the grid is not an axis-aligned shoe-box?");
285 double relevant_coord = fcent[pddir];
286 double pressure = pdrop*(1.0 - (relevant_coord - cmin)/cdelta);
287 int bid = bface_iters[i]->boundaryId();
288 fbcs.setCanonicalBoundaryId(bid, canon_pos + 1);
291 if (twodim_hack && canon_pos >= 4) {
@ Neumann
Definition: BoundaryConditions.hpp:58
@ Dirichlet
Definition: BoundaryConditions.hpp:58
bool isNeumann() const
Type query.
Definition: BoundaryConditions.hpp:91
A class for representing a flow boundary condition.
Definition: BoundaryConditions.hpp:122
double outflux() const
Query a Neumann condition.
Definition: BoundaryConditions.hpp:154
A class for representing a saturation boundary condition.
Definition: BoundaryConditions.hpp:176
Definition: BlackoilFluid.hpp:32
void storeFlowCond(BCs &bcs, const std::vector< BoundaryFaceInfo > &bfinfo, const std::array< FlowBC, 6 > &fconditions, const std::array< double, 6 > &side_areas)
Definition: PeriodicHelpers.hpp:54
void storeSatCond(BCs &bcs, const std::vector< BoundaryFaceInfo > &bfinfo, const std::array< SatBC, 6 > &sconditions)
Definition: PeriodicHelpers.hpp:74
void createPeriodicImpl(BCs &fbcs, std::vector< BoundaryFaceInfo > &bfinfo, std::array< double, 6 > &side_areas, const GridInterface &g, const std::array< bool, 2 *GridInterface::Dimension > &is_periodic, double spatial_tolerance=1e-6)
Common implementation for the various createPeriodic functions.
Definition: PeriodicHelpers.hpp:183
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
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
void findPeriodicPartners(std::vector< BoundaryFaceInfo > &bfinfo, std::array< double, 6 > &side_areas, const GridView &g, const std::array< bool, 2 *GridView::dimension > &is_periodic, double spatial_tolerance=1e-6)
Common implementation for the various createPeriodic functions.
Definition: BoundaryPeriodicity.hpp:86
std::array< bool, 6 > extractPeriodic(const std::array< BC, 6 > &bcs)
Definition: PeriodicHelpers.hpp:86