127 std::cout <<
"==================== Unused parameters: ====================\n";
128 param.displayUsage();
129 std::cout <<
"================================================================\n";
135 template <
class Traits>
137 BoundaryConditionType bctype,
138 double perm_threshold,
139 double residual_tolerance,
140 int linsolver_verbosity,
144 double linsolver_prolongate_factor,
145 int linsolver_smooth_steps,
146 const double gravity)
149 residual_tolerance_ = residual_tolerance;
150 linsolver_verbosity_ = linsolver_verbosity;
151 linsolver_type_ = linsolver_type;
152 linsolver_maxit_ = linsolver_maxit;
153 linsolver_prolongate_factor_ = linsolver_prolongate_factor;
154 linsolver_smooth_steps_ = linsolver_smooth_steps;
155 twodim_hack_ = twodim_hack;
160 bool periodic_ext = (bctype_ == Periodic);
161 bool turn_normals =
false;
162 bool clip_z = (bctype_ == Linear || bctype_ == Periodic);
163 bool unique_bids = (bctype_ == Linear || bctype_ == Periodic);
164 std::string rock_list(
"no_list");
166 periodic_ext, turn_normals, clip_z, unique_bids,
167 perm_threshold, rock_list,
168 useJ<ResProp>(), 1.0, 0.0,
170 ginterf_.init(grid_);
176 template <
class Traits>
177 inline const typename UpscalerBase<Traits>::GridType&
186 template <
class Traits>
190 if ((type == Periodic && bctype_ != Periodic)
191 || (type != Periodic && bctype_ == Periodic)) {
192 OPM_THROW(std::runtime_error,
"Cannot switch to or from Periodic boundary condition, "
193 "periodic must be set in init() params.");
196 if (type == Periodic || type == Linear) {
197 grid_.setUniqueBoundaryIds(
true);
199 grid_.setUniqueBoundaryIds(
false);
207 template <
class Traits>
211 res_prop_.permeabilityModifiable(cell_index) = k;
217 template <
class Traits>
222 return upscaleEffectivePerm(fluid);
228 template <
class Traits>
229 template <
class Flu
idInterface>
233 int num_cells = ginterf_.numberOfCells();
235 std::vector<double> src(num_cells, 0.0);
237 std::vector<double> sat(num_cells, 1.0);
239 Dune::FieldVector<double, 3> gravity(0.0);
240 gravity[2] = gravity_;
242 permtensor_t upscaled_K(3, 3, (
double*)0);
243 for (
int pdd = 0; pdd < Dimension; ++pdd) {
249 flow_solver_.init(ginterf_, res_prop_, gravity, bcond_);
253 bool same_matrix = (bctype_ != Fixed) && (pdd != 0);
254 flow_solver_.solve(fluid, sat, bcond_, src, residual_tolerance_,
255 linsolver_verbosity_,
256 linsolver_type_, same_matrix,
257 linsolver_maxit_, linsolver_prolongate_factor_,
258 linsolver_smooth_steps_);
259 double max_mod = flow_solver_.postProcessFluxes();
260 std::cout <<
"Max mod = " << max_mod << std::endl;
263 double Q[Dimension] = { 0 };
266 Q[pdd] = computeAverageVelocity(flow_solver_.getSolution(), pdd, pdd);
270 for (
int i = 0; i < Dimension; ++i) {
271 Q[i] = computeAverageVelocity(flow_solver_.getSolution(), i, pdd);
275 OPM_THROW(std::runtime_error,
"Unknown boundary type: " + std::to_string(bctype_));
277 double delta = computeDelta(pdd);
278 for (
int i = 0; i < Dimension; ++i) {
279 upscaled_K(i, pdd) = Q[i] * delta;
288 template <
class Traits>
289 template <
class FlowSol>
290 double UpscalerBase<Traits>::computeAverageVelocity(
const FlowSol& flow_solution,
292 const int pdrop_dir)
const
294 double side1_flux = 0.0;
295 double side2_flux = 0.0;
296 double side1_area = 0.0;
297 double side2_area = 0.0;
299 for (CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
300 for (FaceIter f = c->facebegin(); f != c->faceend(); ++f) {
302 int canon_bid = bcond_.getCanonicalBoundaryId(f->boundaryId());
303 if ((canon_bid - 1)/2 == flow_dir) {
304 double flux = flow_solution.outflux(f);
305 double area = f->area();
306 double norm_comp = f->normal()[flow_dir];
308 if (canon_bid - 1 == 2*flow_dir) {
309 if (flow_dir == pdrop_dir && flux > 0.0) {
311 std::cerr <<
"Flow may be in wrong direction at bid: " << f->boundaryId()<<
" (canonical: "<<canon_bid
312 <<
") Magnitude: " << std::fabs(flux) << std::endl;
316 side1_flux += flux*norm_comp;
319 assert(canon_bid - 1 == 2*flow_dir + 1);
320 if (flow_dir == pdrop_dir && flux < 0.0) {
322 std::cerr <<
"Flow may be in wrong direction at bid: " << f->boundaryId()
323 <<
" Magnitude: " << std::fabs(flux) << std::endl;
327 side2_flux += flux*norm_comp;
335 return 0.5*(side1_flux/side1_area + side2_flux/side2_area);
341 template <
class Traits>
342 inline double UpscalerBase<Traits>::computeDelta(
const int flow_dir)
const
344 double side1_pos = 0.0;
345 double side2_pos = 0.0;
346 double side1_area = 0.0;
347 double side2_area = 0.0;
348 for (CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
349 for (FaceIter f = c->facebegin(); f != c->faceend(); ++f) {
351 int canon_bid = bcond_.getCanonicalBoundaryId(f->boundaryId());
352 if ((canon_bid - 1)/2 == flow_dir) {
353 double area = f->area();
354 double pos_comp = f->centroid()[flow_dir];
355 if (canon_bid - 1 == 2*flow_dir) {
356 side1_pos += area*pos_comp;
359 side2_pos += area*pos_comp;
367 return side2_pos/side2_area - side1_pos/side1_area;
373 template <
class Traits>
376 double total_vol = 0.0;
377 double total_pore_vol = 0.0;
378 for (
CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
379 total_vol += c->volume();
380 total_pore_vol += c->volume()*res_prop_.porosity(c->index());
382 return total_pore_vol/total_vol;
386 template <
class Traits>
389 double total_net_vol = 0.0;
390 double total_pore_vol = 0.0;
391 for (
CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
392 total_net_vol += c->volume()*res_prop_.ntg(c->index());
393 total_pore_vol += c->volume()*res_prop_.porosity(c->index())*res_prop_.ntg(c->index());
395 if (total_net_vol>0.0)
return total_pore_vol/total_net_vol;
399 template <
class Traits>
402 double total_vol = 0.0;
403 double total_net_vol = 0.0;
404 for (
CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
405 total_vol += c->volume();
406 total_net_vol += c->volume()*res_prop_.ntg(c->index());
408 return total_net_vol/total_vol;
411 template <
class Traits>
414 double total_swcr = 0.0;
415 double total_pore_vol = 0.0;
417 for (
CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
418 total_swcr += c->volume()*res_prop_.porosity(c->index())*res_prop_.ntg(c->index())*res_prop_.swcr(c->index());
419 total_pore_vol += c->volume()*res_prop_.porosity(c->index())*res_prop_.ntg(c->index());
423 for (
CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
424 total_swcr += c->volume()*res_prop_.porosity(c->index())*res_prop_.swcr(c->index());
425 total_pore_vol += c->volume()*res_prop_.porosity(c->index());
428 return total_swcr/total_pore_vol;
431 template <
class Traits>
434 double total_sowcr = 0.0;
435 double total_pore_vol = 0.0;
437 for (
CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
438 total_sowcr += c->volume()*res_prop_.porosity(c->index())*res_prop_.ntg(c->index())*res_prop_.sowcr(c->index());
439 total_pore_vol += c->volume()*res_prop_.porosity(c->index())*res_prop_.ntg(c->index());
443 for (
CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
444 total_sowcr += c->volume()*res_prop_.porosity(c->index())*res_prop_.sowcr(c->index());
445 total_pore_vol += c->volume()*res_prop_.porosity(c->index());
448 return total_sowcr/total_pore_vol;