35 #ifndef OPM_UPSCALERBASE_IMPL_HEADER
36 #define OPM_UPSCALERBASE_IMPL_HEADER
38 #include <opm/porsol/common/setupGridAndProps.hpp>
39 #include <opm/porsol/common/setupBoundaryConditions.hpp>
40 #include <opm/porsol/common/ReservoirPropertyTracerFluid.hpp>
47 template <
class Traits>
51 residual_tolerance_(1e-8),
53 linsolver_prolongate_factor_(1.0),
54 linsolver_verbosity_(0),
56 linsolver_smooth_steps_(1)
63 template <
class Traits>
73 template <
class Traits>
81 int bct = param.get<
int>(
"boundary_condition_type");
88 twodim_hack_ = param.getDefault(
"2d_hack", twodim_hack_);
89 residual_tolerance_ = param.getDefault(
"residual_tolerance", residual_tolerance_);
90 linsolver_verbosity_ = param.getDefault(
"linsolver_verbosity", linsolver_verbosity_);
91 linsolver_type_ = param.getDefault(
"linsolver_type", linsolver_type_);
92 linsolver_maxit_ = param.getDefault(
"linsolver_max_iterations", linsolver_maxit_);
93 linsolver_prolongate_factor_ = param.getDefault(
"linsolver_prolongate_factor", linsolver_prolongate_factor_);
94 linsolver_smooth_steps_ = param.getDefault(
"linsolver_smooth_steps", linsolver_smooth_steps_);
99 Opm::parameter::ParameterGroup temp_param = param;
100 if (bctype_ == Linear || bctype_ == Periodic) {
101 if (!temp_param.has(
"use_unique_boundary_ids")) {
102 temp_param.insertParameter(
"use_unique_boundary_ids",
"true");
105 if (bctype_ == Periodic) {
106 if (!temp_param.has(
"periodic_extension")) {
107 temp_param.insertParameter(
"periodic_extension",
"true");
111 setupGridAndProps(temp_param, grid_, res_prop_);
112 ginterf_.init(grid_);
118 template <
class Traits>
122 std::cout <<
"==================== Unused parameters: ====================\n";
123 param.displayUsage();
124 std::cout <<
"================================================================\n";
130 template <
class Traits>
133 double perm_threshold,
134 double residual_tolerance,
135 int linsolver_verbosity,
139 double linsolver_prolongate_factor,
140 int linsolver_smooth_steps)
143 residual_tolerance_ = residual_tolerance;
144 linsolver_verbosity_ = linsolver_verbosity;
145 linsolver_type_ = linsolver_type;
146 linsolver_maxit_ = linsolver_maxit;
147 linsolver_prolongate_factor_ = linsolver_prolongate_factor;
148 linsolver_smooth_steps_ = linsolver_smooth_steps;
149 twodim_hack_ = twodim_hack;
152 bool periodic_ext = (bctype_ == Periodic);
153 bool turn_normals =
false;
154 bool clip_z = (bctype_ == Periodic);
155 bool unique_bids = (bctype_ == Linear || bctype_ == Periodic);
156 std::string rock_list(
"no_list");
157 setupGridAndPropsEclipse(deck,
158 periodic_ext, turn_normals, clip_z, unique_bids,
159 perm_threshold, rock_list,
160 useJ<ResProp>(), 1.0, 0.0,
162 ginterf_.init(grid_);
168 template <
class Traits>
178 template <
class Traits>
182 if ((type == Periodic && bctype_ != Periodic)
183 || (type != Periodic && bctype_ == Periodic)) {
184 OPM_THROW(std::runtime_error,
"Cannot switch to or from Periodic boundary condition, "
185 "periodic must be set in init() params.");
188 if (type == Periodic || type == Linear) {
189 grid_.setUniqueBoundaryIds(
true);
191 grid_.setUniqueBoundaryIds(
false);
199 template <
class Traits>
203 res_prop_.permeabilityModifiable(cell_index) = k;
209 template <
class Traits>
213 ReservoirPropertyTracerFluid fluid;
214 return upscaleEffectivePerm(fluid);
220 template <
class Traits>
221 template <
class Flu
idInterface>
225 int num_cells = ginterf_.numberOfCells();
227 std::vector<double> src(num_cells, 0.0);
229 std::vector<double> sat(num_cells, 1.0);
231 Dune::FieldVector<double, 3> gravity(0.0);
235 for (
int pdd = 0; pdd < Dimension; ++pdd) {
236 setupUpscalingConditions(ginterf_, bctype_, pdd, 1.0, 1.0, twodim_hack_, bcond_);
241 flow_solver_.init(ginterf_, res_prop_, gravity, bcond_);
245 bool same_matrix = (bctype_ != Fixed) && (pdd != 0);
246 flow_solver_.solve(fluid, sat, bcond_, src, residual_tolerance_,
247 linsolver_verbosity_,
248 linsolver_type_, same_matrix,
249 linsolver_maxit_, linsolver_prolongate_factor_,
250 linsolver_smooth_steps_);
251 double max_mod = flow_solver_.postProcessFluxes();
252 std::cout <<
"Max mod = " << max_mod << std::endl;
255 double Q[Dimension] = { 0 };
258 Q[pdd] = computeAverageVelocity(flow_solver_.getSolution(), pdd, pdd);
262 for (
int i = 0; i < Dimension; ++i) {
263 Q[i] = computeAverageVelocity(flow_solver_.getSolution(), i, pdd);
267 OPM_THROW(std::runtime_error,
"Unknown boundary type: " << bctype_);
269 double delta = computeDelta(pdd);
270 for (
int i = 0; i < Dimension; ++i) {
271 upscaled_K(i, pdd) = Q[i] * delta;
280 template <
class Traits>
281 template <
class FlowSol>
284 const int pdrop_dir)
const
286 double side1_flux = 0.0;
287 double side2_flux = 0.0;
288 double side1_area = 0.0;
289 double side2_area = 0.0;
292 int num_bdyfaces = 0;
296 for (
CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
297 for (
FaceIter f = c->facebegin(); f != c->faceend(); ++f) {
301 int canon_bid = bcond_.getCanonicalBoundaryId(f->boundaryId());
302 if ((canon_bid - 1)/2 == flow_dir) {
303 double flux = flow_solution.outflux(f);
304 double area = f->area();
305 double norm_comp = f->normal()[flow_dir];
307 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);
321 if (flow_dir == pdrop_dir && flux < 0.0) {
323 std::cerr <<
"Flow may be in wrong direction at bid: " << f->boundaryId()
324 <<
" Magnitude: " << std::fabs(flux) << std::endl;
328 side2_flux += flux*norm_comp;
338 return 0.5*(side1_flux/side1_area + side2_flux/side2_area);
344 template <
class Traits>
347 double side1_pos = 0.0;
348 double side2_pos = 0.0;
349 double side1_area = 0.0;
350 double side2_area = 0.0;
351 for (
CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
352 for (
FaceIter f = c->facebegin(); f != c->faceend(); ++f) {
354 int canon_bid = bcond_.getCanonicalBoundaryId(f->boundaryId());
355 if ((canon_bid - 1)/2 == flow_dir) {
356 double area = f->area();
357 double pos_comp = f->centroid()[flow_dir];
358 if (canon_bid - 1 == 2*flow_dir) {
359 side1_pos += area*pos_comp;
362 side2_pos += area*pos_comp;
370 return side2_pos/side2_area - side1_pos/side1_area;
376 template <
class Traits>
379 double total_vol = 0.0;
380 double total_pore_vol = 0.0;
381 for (
CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
382 total_vol += c->volume();
383 total_pore_vol += c->volume()*res_prop_.porosity(c->index());
385 return total_pore_vol/total_vol;
389 template <
class Traits>
392 double total_net_vol = 0.0;
393 double total_pore_vol = 0.0;
394 for (
CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
395 total_net_vol += c->volume()*res_prop_.ntg(c->index());
396 total_pore_vol += c->volume()*res_prop_.porosity(c->index())*res_prop_.ntg(c->index());
398 if (total_net_vol>0.0)
return total_pore_vol/total_net_vol;
402 template <
class Traits>
405 double total_vol = 0.0;
406 double total_net_vol = 0.0;
407 for (
CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
408 total_vol += c->volume();
409 total_net_vol += c->volume()*res_prop_.ntg(c->index());
411 return total_net_vol/total_vol;
414 template <
class Traits>
417 double total_swcr = 0.0;
418 double total_pore_vol = 0.0;
420 for (
CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
421 total_swcr += c->volume()*res_prop_.porosity(c->index())*res_prop_.ntg(c->index())*res_prop_.swcr(c->index());
422 total_pore_vol += c->volume()*res_prop_.porosity(c->index())*res_prop_.ntg(c->index());
426 for (
CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
427 total_swcr += c->volume()*res_prop_.porosity(c->index())*res_prop_.swcr(c->index());
428 total_pore_vol += c->volume()*res_prop_.porosity(c->index());
431 return total_swcr/total_pore_vol;
434 template <
class Traits>
437 double total_sowcr = 0.0;
438 double total_pore_vol = 0.0;
440 for (
CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
441 total_sowcr += c->volume()*res_prop_.porosity(c->index())*res_prop_.ntg(c->index())*res_prop_.sowcr(c->index());
442 total_pore_vol += c->volume()*res_prop_.porosity(c->index())*res_prop_.ntg(c->index());
446 for (
CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
447 total_sowcr += c->volume()*res_prop_.porosity(c->index())*res_prop_.sowcr(c->index());
448 total_pore_vol += c->volume()*res_prop_.porosity(c->index());
451 return total_sowcr/total_pore_vol;
458 #endif // OPM_UPSCALERBASE_IMPL_HEADER
Dune::CpGrid GridType
Definition: UpscalerBase.hpp:57
Definition: applier.hpp:18
virtual void initFinal(const Opm::parameter::ParameterGroup ¶m)
Definition: UpscalerBase_impl.hpp:119
const GridType & grid() const
Access the grid.
Definition: UpscalerBase_impl.hpp:170
ResProp::MutablePermTensor permtensor_t
A type for the upscaled permeability.
Definition: UpscalerBase.hpp:63
double upscaleNetPorosity() const
Definition: UpscalerBase_impl.hpp:390
void setBoundaryConditionType(BoundaryConditionType type)
Definition: UpscalerBase_impl.hpp:180
double computeDelta(const int flow_dir) const
Definition: UpscalerBase_impl.hpp:345
CellIter::FaceIterator FaceIter
Definition: UpscalerBase.hpp:128
permtensor_t upscaleEffectivePerm(const FluidInterface &fluid)
Definition: UpscalerBase_impl.hpp:223
virtual void initImpl(const Opm::parameter::ParameterGroup ¶m)
Definition: UpscalerBase_impl.hpp:74
double upscaleNTG() const
Definition: UpscalerBase_impl.hpp:403
double upscalePorosity() const
Definition: UpscalerBase_impl.hpp:377
Opm::DeckConstPtr deck(parser->parseFile(file, parseMode))
UpscalerBase()
Default constructor.
Definition: UpscalerBase_impl.hpp:48
double computeAverageVelocity(const FlowSol &flow_solution, const int flow_dir, const int pdrop_dir) const
Definition: UpscalerBase_impl.hpp:282
BoundaryConditionType
Definition: UpscalerBase.hpp:65
void init(const Opm::parameter::ParameterGroup ¶m)
Initializes the upscaler from parameters.
Definition: UpscalerBase_impl.hpp:64
GridInterface::CellIterator CellIter
Definition: UpscalerBase.hpp:127
permtensor_t upscaleSinglePhase()
Definition: UpscalerBase_impl.hpp:211
double upscaleSWCR(const bool NTG) const
Definition: UpscalerBase_impl.hpp:415
double upscaleSOWCR(const bool NTG) const
Definition: UpscalerBase_impl.hpp:435
void setPermeability(const int cell_index, const permtensor_t &k)
Definition: UpscalerBase_impl.hpp:201