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");
82 bctype_ =
static_cast<BoundaryConditionType
>(bct);
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::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");
104 if (!temp_param.has(
"clip_z")) {
105 temp_param.insertParameter(
"clip_z",
"true");
108 if (bctype_ == Periodic) {
109 if (!temp_param.has(
"periodic_extension")) {
110 temp_param.insertParameter(
"periodic_extension",
"true");
115 ginterf_.init(grid_);
121 template <
class Traits>
122 inline void UpscalerBase<Traits>::initFinal(
const Opm::ParameterGroup& param)
125 std::cout <<
"==================== Unused parameters: ====================\n";
126 param.displayUsage();
127 std::cout <<
"================================================================\n";
133 template <
class Traits>
135 BoundaryConditionType bctype,
136 double perm_threshold,
137 double residual_tolerance,
138 int linsolver_verbosity,
142 double linsolver_prolongate_factor,
143 int linsolver_smooth_steps,
144 const double gravity)
147 residual_tolerance_ = residual_tolerance;
148 linsolver_verbosity_ = linsolver_verbosity;
149 linsolver_type_ = linsolver_type;
150 linsolver_maxit_ = linsolver_maxit;
151 linsolver_prolongate_factor_ = linsolver_prolongate_factor;
152 linsolver_smooth_steps_ = linsolver_smooth_steps;
153 twodim_hack_ = twodim_hack;
158 bool periodic_ext = (bctype_ == Periodic);
159 bool turn_normals =
false;
160 bool clip_z = (bctype_ == Linear || bctype_ == Periodic);
161 bool unique_bids = (bctype_ == Linear || bctype_ == Periodic);
162 std::string rock_list(
"no_list");
164 periodic_ext, turn_normals, clip_z, unique_bids,
165 perm_threshold, rock_list,
166 useJ<ResProp>(), 1.0, 0.0,
168 ginterf_.init(grid_);
174 template <
class Traits>
175 inline const typename UpscalerBase<Traits>::GridType&
184 template <
class Traits>
188 if ((type == Periodic && bctype_ != Periodic)
189 || (type != Periodic && bctype_ == Periodic)) {
190 OPM_THROW(std::runtime_error,
"Cannot switch to or from Periodic boundary condition, "
191 "periodic must be set in init() params.");
194 if (type == Periodic || type == Linear) {
195 grid_.setUniqueBoundaryIds(
true);
197 grid_.setUniqueBoundaryIds(
false);
205 template <
class Traits>
209 res_prop_.permeabilityModifiable(cell_index) = k;
215 template <
class Traits>
220 return upscaleEffectivePerm(fluid);
226 template <
class Traits>
227 template <
class Flu
idInterface>
231 int num_cells = ginterf_.numberOfCells();
233 std::vector<double> src(num_cells, 0.0);
235 std::vector<double> sat(num_cells, 1.0);
237 Dune::FieldVector<double, 3> gravity(0.0);
238 gravity[2] = gravity_;
240 permtensor_t upscaled_K(3, 3, (
double*)0);
241 for (
int pdd = 0; pdd < Dimension; ++pdd) {
247 flow_solver_.init(ginterf_, res_prop_, gravity, bcond_);
251 bool same_matrix = (bctype_ != Fixed) && (pdd != 0);
252 flow_solver_.solve(fluid, sat, bcond_, src, residual_tolerance_,
253 linsolver_verbosity_,
254 linsolver_type_, same_matrix,
255 linsolver_maxit_, linsolver_prolongate_factor_,
256 linsolver_smooth_steps_);
257 double max_mod = flow_solver_.postProcessFluxes();
258 std::cout <<
"Max mod = " << max_mod << std::endl;
261 double Q[Dimension] = { 0 };
264 Q[pdd] = computeAverageVelocity(flow_solver_.getSolution(), pdd, pdd);
268 for (
int i = 0; i < Dimension; ++i) {
269 Q[i] = computeAverageVelocity(flow_solver_.getSolution(), i, pdd);
273 OPM_THROW(std::runtime_error,
"Unknown boundary type: " << bctype_);
275 double delta = computeDelta(pdd);
276 for (
int i = 0; i < Dimension; ++i) {
277 upscaled_K(i, pdd) = Q[i] * delta;
286 template <
class Traits>
287 template <
class FlowSol>
288 double UpscalerBase<Traits>::computeAverageVelocity(
const FlowSol& flow_solution,
290 const int pdrop_dir)
const
292 double side1_flux = 0.0;
293 double side2_flux = 0.0;
294 double side1_area = 0.0;
295 double side2_area = 0.0;
298 int num_bdyfaces = 0;
302 for (CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
303 for (FaceIter f = c->facebegin(); f != c->faceend(); ++f) {
307 int canon_bid = bcond_.getCanonicalBoundaryId(f->boundaryId());
308 if ((canon_bid - 1)/2 == flow_dir) {
309 double flux = flow_solution.outflux(f);
310 double area = f->area();
311 double norm_comp = f->normal()[flow_dir];
313 if (canon_bid - 1 == 2*flow_dir) {
315 if (flow_dir == pdrop_dir && flux > 0.0) {
317 std::cerr <<
"Flow may be in wrong direction at bid: " << f->boundaryId()<<
" (canonical: "<<canon_bid
318 <<
") Magnitude: " << std::fabs(flux) << std::endl;
322 side1_flux += flux*norm_comp;
325 assert(canon_bid - 1 == 2*flow_dir + 1);
327 if (flow_dir == pdrop_dir && flux < 0.0) {
329 std::cerr <<
"Flow may be in wrong direction at bid: " << f->boundaryId()
330 <<
" Magnitude: " << std::fabs(flux) << std::endl;
334 side2_flux += flux*norm_comp;
344 return 0.5*(side1_flux/side1_area + side2_flux/side2_area);
350 template <
class Traits>
351 inline double UpscalerBase<Traits>::computeDelta(
const int flow_dir)
const
353 double side1_pos = 0.0;
354 double side2_pos = 0.0;
355 double side1_area = 0.0;
356 double side2_area = 0.0;
357 for (CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
358 for (FaceIter f = c->facebegin(); f != c->faceend(); ++f) {
360 int canon_bid = bcond_.getCanonicalBoundaryId(f->boundaryId());
361 if ((canon_bid - 1)/2 == flow_dir) {
362 double area = f->area();
363 double pos_comp = f->centroid()[flow_dir];
364 if (canon_bid - 1 == 2*flow_dir) {
365 side1_pos += area*pos_comp;
368 side2_pos += area*pos_comp;
376 return side2_pos/side2_area - side1_pos/side1_area;
382 template <
class Traits>
385 double total_vol = 0.0;
386 double total_pore_vol = 0.0;
387 for (
CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
388 total_vol += c->volume();
389 total_pore_vol += c->volume()*res_prop_.porosity(c->index());
391 return total_pore_vol/total_vol;
395 template <
class Traits>
398 double total_net_vol = 0.0;
399 double total_pore_vol = 0.0;
400 for (
CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
401 total_net_vol += c->volume()*res_prop_.ntg(c->index());
402 total_pore_vol += c->volume()*res_prop_.porosity(c->index())*res_prop_.ntg(c->index());
404 if (total_net_vol>0.0)
return total_pore_vol/total_net_vol;
408 template <
class Traits>
411 double total_vol = 0.0;
412 double total_net_vol = 0.0;
413 for (
CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
414 total_vol += c->volume();
415 total_net_vol += c->volume()*res_prop_.ntg(c->index());
417 return total_net_vol/total_vol;
420 template <
class Traits>
423 double total_swcr = 0.0;
424 double total_pore_vol = 0.0;
426 for (
CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
427 total_swcr += c->volume()*res_prop_.porosity(c->index())*res_prop_.ntg(c->index())*res_prop_.swcr(c->index());
428 total_pore_vol += c->volume()*res_prop_.porosity(c->index())*res_prop_.ntg(c->index());
432 for (
CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
433 total_swcr += c->volume()*res_prop_.porosity(c->index())*res_prop_.swcr(c->index());
434 total_pore_vol += c->volume()*res_prop_.porosity(c->index());
437 return total_swcr/total_pore_vol;
440 template <
class Traits>
443 double total_sowcr = 0.0;
444 double total_pore_vol = 0.0;
446 for (
CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
447 total_sowcr += c->volume()*res_prop_.porosity(c->index())*res_prop_.ntg(c->index())*res_prop_.sowcr(c->index());
448 total_pore_vol += c->volume()*res_prop_.porosity(c->index())*res_prop_.ntg(c->index());
452 for (
CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
453 total_sowcr += c->volume()*res_prop_.porosity(c->index())*res_prop_.sowcr(c->index());
454 total_pore_vol += c->volume()*res_prop_.porosity(c->index());
457 return total_sowcr/total_pore_vol;
464 #endif // OPM_UPSCALERBASE_IMPL_HEADER
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:151
void setupUpscalingConditions(const GridInterface &g, int bct, int pddir, double pdrop, double bdy_sat, bool twodim_hack, BCs &bcs)
Definition: setupBoundaryConditions.hpp:99
Definition: GridInterfaceEuler.hpp:376
A base class for upscaling.
Definition: UpscalerBase.hpp:55
void setPermeability(const int cell_index, const permtensor_t &k)
Set the permeability of a cell directly.
Definition: UpscalerBase_impl.hpp:207
void setBoundaryConditionType(BoundaryConditionType type)
Set boundary condition type.
Definition: UpscalerBase_impl.hpp:186
const GridType & grid() const
Access the grid.
Definition: UpscalerBase_impl.hpp:176
void setupGridAndProps(const Opm::ParameterGroup ¶m, Dune::CpGrid &grid, ResProp< 3 > &res_prop)
Definition: setupGridAndProps.hpp:70
double upscaleNetPorosity() const
Compute upscaled net porosity.
Definition: UpscalerBase_impl.hpp:396
Definition: ReservoirPropertyTracerFluid.hpp:40
double upscalePorosity() const
Compute upscaled porosity.
Definition: UpscalerBase_impl.hpp:383
UpscalerBase()
Default constructor.
Definition: UpscalerBase_impl.hpp:48
double upscaleSOWCR(const bool NTG) const
Compute upscaled SOWCR.
Definition: UpscalerBase_impl.hpp:441
ResProp::MutablePermTensor permtensor_t
A type for the upscaled permeability.
Definition: UpscalerBase.hpp:66
double upscaleNTG() const
Compute upscaled NTG.
Definition: UpscalerBase_impl.hpp:409
void init(const Opm::ParameterGroup ¶m)
Initializes the upscaler from parameters.
Definition: UpscalerBase_impl.hpp:64
permtensor_t upscaleSinglePhase()
Does a single-phase upscaling.
Definition: UpscalerBase_impl.hpp:217
double upscaleSWCR(const bool NTG) const
Compute upscaled SWCR.
Definition: UpscalerBase_impl.hpp:421