35#ifndef OPM_UPSCALERBASE_IMPL_HEADER
36#define OPM_UPSCALERBASE_IMPL_HEADER
47 template <
class Traits>
51 residual_tolerance_(1e-8),
53 linsolver_prolongate_factor_(1.0),
54 linsolver_verbosity_(0),
56 linsolver_smooth_steps_(1),
64 template <
class Traits>
74 template <
class Traits>
82 int bct = param.get<
int>(
"boundary_condition_type");
89 twodim_hack_ = param.getDefault(
"2d_hack", twodim_hack_);
90 residual_tolerance_ = param.getDefault(
"residual_tolerance", residual_tolerance_);
91 linsolver_verbosity_ = param.getDefault(
"linsolver_verbosity", linsolver_verbosity_);
92 linsolver_type_ = param.getDefault(
"linsolver_type", linsolver_type_);
93 linsolver_maxit_ = param.getDefault(
"linsolver_max_iterations", linsolver_maxit_);
94 linsolver_prolongate_factor_ = param.getDefault(
"linsolver_prolongate_factor", linsolver_prolongate_factor_);
95 linsolver_smooth_steps_ = param.getDefault(
"linsolver_smooth_steps", linsolver_smooth_steps_);
96 gravity_ = param.getDefault(
"gravity", gravity_);
101 Opm::ParameterGroup temp_param = param;
102 if (bctype_ == Linear || bctype_ == Periodic) {
103 if (!temp_param.has(
"use_unique_boundary_ids")) {
104 temp_param.insertParameter(
"use_unique_boundary_ids",
"true");
106 if (!temp_param.has(
"clip_z")) {
107 temp_param.insertParameter(
"clip_z",
"true");
110 if (bctype_ == Periodic) {
111 if (!temp_param.has(
"periodic_extension")) {
112 temp_param.insertParameter(
"periodic_extension",
"true");
117 ginterf_.init(grid_);
123 template <
class Traits>
127 std::cout <<
"==================== Unused parameters: ====================\n";
128 param.displayUsage();
129 std::cout <<
"================================================================\n";
135 template <
class Traits>
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>
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_;
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>
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>
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;
Definition: GridInterfaceEuler.hpp:368
Intersection (face) iterator for solver-near grid interface.
Definition: GridInterfaceEuler.hpp:241
Definition: ReservoirPropertyTracerFluid.hpp:40
double computeDelta(const int flow_dir) const
Definition: UpscalerBase_impl.hpp:342
permtensor_t upscaleSinglePhase()
Definition: UpscalerBase_impl.hpp:219
void setPermeability(const int cell_index, const permtensor_t &k)
Definition: UpscalerBase_impl.hpp:209
double computeAverageVelocity(const FlowSol &flow_solution, const int flow_dir, const int pdrop_dir) const
Definition: UpscalerBase_impl.hpp:290
double upscalePorosity() const
Definition: UpscalerBase_impl.hpp:374
ResProp::MutablePermTensor permtensor_t
A type for the upscaled permeability.
Definition: UpscalerBase.hpp:66
const GridType & grid() const
Access the grid.
Definition: UpscalerBase_impl.hpp:178
double upscaleNTG() const
Definition: UpscalerBase_impl.hpp:400
virtual void initFinal(const Opm::ParameterGroup ¶m)
Definition: UpscalerBase_impl.hpp:124
BoundaryConditionType
Definition: UpscalerBase.hpp:68
permtensor_t upscaleEffectivePerm(const FluidInterface &fluid)
Definition: UpscalerBase_impl.hpp:231
Dune::CpGrid GridType
Definition: UpscalerBase.hpp:60
double upscaleSOWCR(const bool NTG) const
Definition: UpscalerBase_impl.hpp:432
virtual void initImpl(const Opm::ParameterGroup ¶m)
Definition: UpscalerBase_impl.hpp:75
void init(const Opm::ParameterGroup ¶m)
Initializes the upscaler from parameters.
Definition: UpscalerBase_impl.hpp:65
void setBoundaryConditionType(BoundaryConditionType type)
Definition: UpscalerBase_impl.hpp:188
double upscaleNetPorosity() const
Definition: UpscalerBase_impl.hpp:387
double upscaleSWCR(const bool NTG) const
Definition: UpscalerBase_impl.hpp:412
UpscalerBase()
Default constructor.
Definition: UpscalerBase_impl.hpp:48
auto deck
Definition: elasticity_upscale_impl.hpp:590
Definition: ImplicitAssembly.hpp:43
void setupGridAndPropsEclipse(const Opm::Deck &deck, bool periodic_extension, bool turn_normals, bool clip_z, bool unique_bids, double perm_threshold, const std::string &rock_list, bool use_jfunction_scaling, double sigma, double theta, Dune::CpGrid &grid, ResProp< 3 > &res_prop)
Definition: setupGridAndProps.hpp:152
void setupUpscalingConditions(const GridInterface &g, int bct, int pddir, double pdrop, double bdy_sat, bool twodim_hack, BCs &bcs)
Definition: setupBoundaryConditions.hpp:99
void setupGridAndProps(const Opm::ParameterGroup ¶m, Dune::CpGrid &grid, ResProp< 3 > &res_prop)
Definition: setupGridAndProps.hpp:71