My Project
PeriodicHelpers.hpp
1 //===========================================================================
2 //
3 // File: PeriodicHelpers.hpp
4 //
5 // Created: Fri Aug 14 10:59:57 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_PERIODICHELPERS_HEADER
37 #define OPENRS_PERIODICHELPERS_HEADER
38 
39 #include "BoundaryPeriodicity.hpp"
40 #include "BoundaryConditions.hpp"
41 #include <dune/common/fvector.hh>
42 #include <array>
43 #include <algorithm>
44 #include <iostream>
45 
46 namespace Opm
47 {
48 
49 
50 
51 
52 
53  template <class BCs>
54  void storeFlowCond(BCs& bcs,
55  const std::vector<BoundaryFaceInfo>& bfinfo,
56  const std::array<FlowBC, 6>& fconditions,
57  const std::array<double, 6>& side_areas)
58  {
59  int num_bdy = bfinfo.size();
60  for (int i = 0; i < num_bdy; ++i) {
61  FlowBC face_bc = fconditions[bfinfo[i].canon_pos];
62  if (face_bc.isNeumann()) {
63  face_bc = FlowBC(FlowBC::Neumann, face_bc.outflux()*bfinfo[i].area/side_areas[bfinfo[i].canon_pos]);
64  }
65  bcs.flowCond(bfinfo[i].bid) = face_bc;
66  }
67  }
68 
69 
70 
71 
72 
73  template <class BCs>
74  void storeSatCond(BCs& bcs,
75  const std::vector<BoundaryFaceInfo>& bfinfo,
76  const std::array<SatBC, 6>& sconditions)
77  {
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];
81  }
82  }
83 
84 
85  template <class BC>
86  std::array<bool, 6> extractPeriodic(const std::array<BC, 6>& bcs)
87  {
88  std::array<bool, 6> retval = {{ bcs[0].isPeriodic(),
89  bcs[1].isPeriodic(),
90  bcs[2].isPeriodic(),
91  bcs[3].isPeriodic(),
92  bcs[4].isPeriodic(),
93  bcs[5].isPeriodic() }};
94  return retval;
95  }
96 
97 
104  template <class BCs, class GridInterface>
105  void createPeriodic(BCs& fbcs,
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)
110  {
111  static_assert(BCs::HasFlowConds, "");
112  static_assert(BCs::HasSatConds, "");
113  // Check the conditions given.
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);
122  }
123  }
124  std::vector<BoundaryFaceInfo> bfinfo;
125  std::array<double, 6> side_areas;
126  createPeriodicImpl(fbcs, bfinfo, side_areas, g, extractPeriodic(fconditions), spatial_tolerance);
127  storeFlowCond(fbcs, bfinfo, fconditions, side_areas);
128  storeSatCond(fbcs, bfinfo, sconditions);
129  }
130 
131 
132 
133 
134  template <class BCs, class GridInterface>
135  void createPeriodic(BCs& fbcs,
136  const GridInterface& g,
137  const std::array<FlowBC, 2*GridInterface::Dimension>& fconditions,
138  double spatial_tolerance = 1e-6)
139  {
140  static_assert(BCs::HasFlowConds, "");
141  static_assert(!BCs::HasSatConds, "");
142  // Check the conditions given.
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());
147  }
148  }
149  std::vector<BoundaryFaceInfo> bfinfo;
150  std::array<double, 6> side_areas;
151  createPeriodicImpl(fbcs, bfinfo, side_areas, g, extractPeriodic(fconditions), spatial_tolerance);
152  storeFlowCond(fbcs, bfinfo, fconditions, side_areas);
153  }
154 
155 
156 
157 
158  template <class BCs, class GridInterface>
159  void createPeriodic(BCs& fbcs,
160  const GridInterface& g,
161  const std::array<SatBC, 2*GridInterface::Dimension>& sconditions,
162  double spatial_tolerance = 1e-6)
163  {
164  static_assert(!BCs::HasFlowConds, "");
165  static_assert(BCs::HasSatConds, "");
166  // Check the conditions given.
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());
171  }
172  }
173  std::vector<BoundaryFaceInfo> bfinfo;
174  std::array<double, 6> side_areas;
175  createPeriodicImpl(fbcs, bfinfo, side_areas, g, extractPeriodic(sconditions), spatial_tolerance);
176  storeSatCond(fbcs, bfinfo, sconditions);
177  }
178 
179 
180 
182  template <class BCs, class GridInterface>
183  void createPeriodicImpl(BCs& fbcs,
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)
189  {
190  findPeriodicPartners(bfinfo, side_areas, g.grid().leafGridView(), is_periodic, spatial_tolerance);
191  int num_bdy = bfinfo.size();
192  // This will likely change with boundarySegmentIndex() instead of boundaryId():
193  int max_bid = num_bdy;
194 
195  // Resize the conditions object. We clear it first to make sure it's all defaults.
196  fbcs.clear();
197  fbcs.resize(max_bid + 1);
198 
199  // Now we have all the info, and the object to store it in. So we store it...
200  for (int i = 0; i < num_bdy; ++i) {
201  int bid1 = bfinfo[i].bid;
202  int bid2 = bfinfo[i].partner_bid;
203  if (bid1 < bid2) {
204  // If there is no periodic partner, bid2 will be zero, so we will not end up here.
205  fbcs.setPeriodicPartners(bid1, bid2);
206  }
207  fbcs.setCanonicalBoundaryId(bid1, bfinfo[i].canon_pos + 1);
208  }
209  }
210 
211 
212 
219  template <class BCs, class GridInterface>
220  void createLinear(BCs& fbcs,
221  const GridInterface& g,
222  const double pdrop,
223  const int pddir,
224  const double bdy_sat,
225  const bool twodim_hack = false,
226  const double spatial_tolerance = 1e-6)
227  {
228  // NOTE: Section copied from createPeriodicImpl().
229  // Pick out all boundary faces, simultaneously find the
230  // bounding box of their centroids, and the max id.
231  typedef typename GridInterface::CellIterator CI;
232  typedef typename CI::FaceIterator FI;
233  typedef typename GridInterface::Vector Vector;
234  std::vector<FI> bface_iters;
235  Vector low(1e100);
236  Vector hi(-1e100);
237  int max_bid = 0;
238  for (CI c = g.cellbegin(); c != g.cellend(); ++c) {
239  for (FI f = c->facebegin(); f != c->faceend(); ++f) {
240  if (f->boundaryId()) {
241  bface_iters.push_back(f);
242  Vector fcent = f->centroid();
243  for (int dd = 0; dd < GridInterface::Dimension; ++dd) {
244  low[dd] = std::min(low[dd], fcent[dd]);
245  hi[dd] = std::max(hi[dd], fcent[dd]);
246  }
247  max_bid = std::max(max_bid, f->boundaryId());
248  }
249  }
250  }
251  int num_bdy = bface_iters.size();
252  if (max_bid != num_bdy) {
253  OPM_THROW(std::runtime_error, "createLinear() assumes that every boundary face has a unique boundary id. That seems to be violated.");
254  }
255  fbcs.resize(max_bid + 1);
256 
257  // Iterate over boundary faces, setting boundary conditions for their corresponding ids.
258  double cmin = low[pddir];
259  double cmax = hi[pddir];
260  double cdelta = cmax - cmin;
261 
262  for (int i = 0; i < num_bdy; ++i) {
263  Vector fcent = bface_iters[i]->centroid();
264  int canon_pos = -1;
265  for (int dd = 0; dd < GridInterface::Dimension; ++dd) {
266  double coord = fcent[dd];
267  if (fabs(coord - low[dd]) <= spatial_tolerance) {
268  canon_pos = 2*dd;
269  break;
270  } else if (fabs(coord - hi[dd]) <= spatial_tolerance) {
271  canon_pos = 2*dd + 1;
272  break;
273  }
274  }
275  if (canon_pos == -1) {
276  std::cerr << "Centroid: " << fcent << "\n";
277  std::cerr << "Bounding box min: " << low << "\n";
278  std::cerr << "Bounding box max: " << hi << "\n";
279  OPM_THROW(std::runtime_error, "Boundary face centroid not on bounding box. Maybe the grid is not an axis-aligned shoe-box?");
280  }
281  double relevant_coord = fcent[pddir];
282  double pressure = pdrop*(1.0 - (relevant_coord - cmin)/cdelta);
283  int bid = bface_iters[i]->boundaryId();
284  fbcs.setCanonicalBoundaryId(bid, canon_pos + 1);
285  fbcs.satCond(bid) = SatBC(SatBC::Dirichlet, bdy_sat);
286  fbcs.flowCond(bid) = FlowBC(FlowBC::Dirichlet, pressure);
287  if (twodim_hack && canon_pos >= 4) {
288  // Noflow for Z- and Z+ boundaries in the 3d-grid-that-is-really-2d case.
289  fbcs.flowCond(bid) = FlowBC(FlowBC::Neumann, 0.0);
290  }
291  }
292  }
293 
294 
295 } // namespace Opm
296 
297 #endif // OPENRS_PERIODICHELPERS_HEADER
A class for representing a flow boundary condition.
Definition: BoundaryConditions.hpp:122
A class for representing a saturation boundary condition.
Definition: BoundaryConditions.hpp:176
Inverting small matrices.
Definition: ImplicitAssembly.hpp:43
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:220
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