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),
64 template <
class Traits>
74 template <
class Traits>
82 int bct = param.get<
int>(
"boundary_condition_type");
83 bctype_ =
static_cast<BoundaryConditionType
>(bct);
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>
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: " << 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;
300 int num_bdyfaces = 0;
304 for (CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
305 for (FaceIter f = c->facebegin(); f != c->faceend(); ++f) {
309 int canon_bid = bcond_.getCanonicalBoundaryId(f->boundaryId());
310 if ((canon_bid - 1)/2 == flow_dir) {
311 double flux = flow_solution.outflux(f);
312 double area = f->area();
313 double norm_comp = f->normal()[flow_dir];
315 if (canon_bid - 1 == 2*flow_dir) {
317 if (flow_dir == pdrop_dir && flux > 0.0) {
319 std::cerr <<
"Flow may be in wrong direction at bid: " << f->boundaryId()<<
" (canonical: "<<canon_bid
320 <<
") Magnitude: " << std::fabs(flux) << std::endl;
324 side1_flux += flux*norm_comp;
327 assert(canon_bid - 1 == 2*flow_dir + 1);
329 if (flow_dir == pdrop_dir && flux < 0.0) {
331 std::cerr <<
"Flow may be in wrong direction at bid: " << f->boundaryId()
332 <<
" Magnitude: " << std::fabs(flux) << std::endl;
336 side2_flux += flux*norm_comp;
346 return 0.5*(side1_flux/side1_area + side2_flux/side2_area);
352 template <
class Traits>
353 inline double UpscalerBase<Traits>::computeDelta(
const int flow_dir)
const
355 double side1_pos = 0.0;
356 double side2_pos = 0.0;
357 double side1_area = 0.0;
358 double side2_area = 0.0;
359 for (CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
360 for (FaceIter f = c->facebegin(); f != c->faceend(); ++f) {
362 int canon_bid = bcond_.getCanonicalBoundaryId(f->boundaryId());
363 if ((canon_bid - 1)/2 == flow_dir) {
364 double area = f->area();
365 double pos_comp = f->centroid()[flow_dir];
366 if (canon_bid - 1 == 2*flow_dir) {
367 side1_pos += area*pos_comp;
370 side2_pos += area*pos_comp;
378 return side2_pos/side2_area - side1_pos/side1_area;
384 template <
class Traits>
387 double total_vol = 0.0;
388 double total_pore_vol = 0.0;
389 for (
CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
390 total_vol += c->volume();
391 total_pore_vol += c->volume()*res_prop_.porosity(c->index());
393 return total_pore_vol/total_vol;
397 template <
class Traits>
400 double total_net_vol = 0.0;
401 double total_pore_vol = 0.0;
402 for (
CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
403 total_net_vol += c->volume()*res_prop_.ntg(c->index());
404 total_pore_vol += c->volume()*res_prop_.porosity(c->index())*res_prop_.ntg(c->index());
406 if (total_net_vol>0.0)
return total_pore_vol/total_net_vol;
410 template <
class Traits>
413 double total_vol = 0.0;
414 double total_net_vol = 0.0;
415 for (
CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
416 total_vol += c->volume();
417 total_net_vol += c->volume()*res_prop_.ntg(c->index());
419 return total_net_vol/total_vol;
422 template <
class Traits>
425 double total_swcr = 0.0;
426 double total_pore_vol = 0.0;
428 for (
CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
429 total_swcr += c->volume()*res_prop_.porosity(c->index())*res_prop_.ntg(c->index())*res_prop_.swcr(c->index());
430 total_pore_vol += c->volume()*res_prop_.porosity(c->index())*res_prop_.ntg(c->index());
434 for (
CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
435 total_swcr += c->volume()*res_prop_.porosity(c->index())*res_prop_.swcr(c->index());
436 total_pore_vol += c->volume()*res_prop_.porosity(c->index());
439 return total_swcr/total_pore_vol;
442 template <
class Traits>
445 double total_sowcr = 0.0;
446 double total_pore_vol = 0.0;
448 for (
CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
449 total_sowcr += c->volume()*res_prop_.porosity(c->index())*res_prop_.ntg(c->index())*res_prop_.sowcr(c->index());
450 total_pore_vol += c->volume()*res_prop_.porosity(c->index())*res_prop_.ntg(c->index());
454 for (
CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
455 total_sowcr += c->volume()*res_prop_.porosity(c->index())*res_prop_.sowcr(c->index());
456 total_pore_vol += c->volume()*res_prop_.porosity(c->index());
459 return total_sowcr/total_pore_vol;
Definition: GridInterfaceEuler.hpp:367
Definition: ReservoirPropertyTracerFluid.hpp:40
A base class for upscaling.
Definition: UpscalerBase.hpp:56
permtensor_t upscaleSinglePhase()
Does a single-phase upscaling.
Definition: UpscalerBase_impl.hpp:219
void setPermeability(const int cell_index, const permtensor_t &k)
Set the permeability of a cell directly.
Definition: UpscalerBase_impl.hpp:209
double upscalePorosity() const
Compute upscaled porosity.
Definition: UpscalerBase_impl.hpp:385
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
Compute upscaled NTG.
Definition: UpscalerBase_impl.hpp:411
double upscaleSOWCR(const bool NTG) const
Compute upscaled SOWCR.
Definition: UpscalerBase_impl.hpp:443
void init(const Opm::ParameterGroup ¶m)
Initializes the upscaler from parameters.
Definition: UpscalerBase_impl.hpp:65
void setBoundaryConditionType(BoundaryConditionType type)
Set boundary condition type.
Definition: UpscalerBase_impl.hpp:188
double upscaleNetPorosity() const
Compute upscaled net porosity.
Definition: UpscalerBase_impl.hpp:398
double upscaleSWCR(const bool NTG) const
Compute upscaled SWCR.
Definition: UpscalerBase_impl.hpp:423
UpscalerBase()
Default constructor.
Definition: UpscalerBase_impl.hpp:48
Inverting small matrices.
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:149
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:69