BlackoilInitialization.hpp
Go to the documentation of this file.
1/*
2 Copyright 2011 SINTEF ICT, Applied Mathematics.
3
4 This file is part of the Open Porous Media project (OPM).
5
6 OPM is free software: you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation, either version 3 of the License, or
9 (at your option) any later version.
10
11 OPM is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with OPM. If not, see <http://www.gnu.org/licenses/>.
18*/
19
20#ifndef OPM_BLACKOILINITIALIZATION_HEADER_INCLUDED
21#define OPM_BLACKOILINITIALIZATION_HEADER_INCLUDED
22
23#include <opm/core/utility/parameters/ParameterGroup.hpp>
24#include <opm/common/ErrorMacros.hpp>
25#include <opm/core/utility/Units.hpp>
26
27#include <iostream>
28
29namespace Opm
30{
31
33 template <class Simulator>
34 class BlackoilInitialization
35 {
36 public:
37 typedef typename Simulator::State State;
38 typedef typename Simulator::Grid Grid;
39 typedef typename Simulator::Fluid Fluid;
40
41 virtual void init(const Opm::parameter::ParameterGroup& param,
42 const Grid& grid,
43 const Fluid& fluid,
44 typename Grid::Vector gravity,
45 State& simstate) = 0;
46 };
47
48
49
51 template <class Simulator>
52 class BasicInitialization : public BlackoilInitialization<Simulator>
53 {
54 public:
55 typedef typename Simulator::State State;
56 typedef typename Simulator::Grid Grid;
57 typedef typename Simulator::Fluid Fluid;
58
59 virtual void init(const Opm::parameter::ParameterGroup& param,
60 const Grid& grid,
61 const Fluid& fluid,
62 typename Grid::Vector gravity,
63 State& simstate)
64 {
65 typedef typename Fluid::CompVec CompVec;
66 typedef typename Fluid::PhaseVec PhaseVec;
67
68 if (param.getDefault("heterogenous_initial_mix", false)) {
69 CompVec init_oil(0.0);
70 init_oil[Fluid::Oil] = 1.0;
71 CompVec init_water(0.0);
72 init_water[Fluid::Water] = 1.0;
73 simstate.cell_z_.resize(grid.numCells());
74 std::fill(simstate.cell_z_.begin(), simstate.cell_z_.begin() + simstate.cell_z_.size()/2, init_oil);
75 std::fill(simstate.cell_z_.begin() + simstate.cell_z_.size()/2, simstate.cell_z_.end(), init_water);
76 OPM_MESSAGE("******* Assuming zero capillary pressures *******");
77 PhaseVec init_p(100.0*Opm::unit::barsa);
78 simstate.cell_pressure_.resize(grid.numCells(), init_p);
79 // if (gravity.two_norm() != 0.0) {
80 // double ref_gravpot = grid.cellCentroid(0)*gravity;
81 // double rho = init_z*fluid_.surfaceDensities(); // Assuming incompressible, and constant initial z.
82 // for (int cell = 1; cell < grid.numCells(); ++cell) {
83 // double press = rho*(grid.cellCentroid(cell)*gravity - ref_gravpot) + simstate.cell_pressure_[0][0];
84 // simstate.cell_pressure_[cell] = PhaseVec(press);
85 // }
86 // }
87 } else if (param.getDefault("unstable_initial_mix", false)) {
88 CompVec init_oil(0.0);
89 init_oil[Fluid::Oil] = 1.0;
90 init_oil[Fluid::Gas] = 0.0;
91 CompVec init_water(0.0);
92 init_water[Fluid::Water] = 1.0;
93 CompVec init_gas(0.0);
94 init_gas[Fluid::Gas] = 150.0;
95 simstate.cell_z_.resize(grid.numCells());
96 std::fill(simstate.cell_z_.begin(),
97 simstate.cell_z_.begin() + simstate.cell_z_.size()/3,
98 init_water);
99 std::fill(simstate.cell_z_.begin() + simstate.cell_z_.size()/3,
100 simstate.cell_z_.begin() + 2*(simstate.cell_z_.size()/3),
101 init_oil);
102 std::fill(simstate.cell_z_.begin() + 2*(simstate.cell_z_.size()/3),
103 simstate.cell_z_.end(),
104 init_gas);
105 OPM_MESSAGE("******* Assuming zero capillary pressures *******");
106 PhaseVec init_p(100.0*Opm::unit::barsa);
107 simstate.cell_pressure_.resize(grid.numCells(), init_p);
108
109 if (gravity.two_norm() != 0.0) {
110
111 typename Fluid::FluidState state = fluid.computeState(simstate.cell_pressure_[0], simstate.cell_z_[0]);
112 simstate.cell_z_[0] *= 1.0/state.total_phase_volume_density_;
113 for (int cell = 1; cell < grid.numCells(); ++cell) {
114 double fluid_vol_dens;
115 int cnt =0;
116 do {
117 double rho = 0.5*((simstate.cell_z_[cell]+simstate.cell_z_[cell-1])*fluid.surfaceDensities());
118 double press = rho*((grid.cellCentroid(cell) - grid.cellCentroid(cell-1))*gravity) + simstate.cell_pressure_[cell-1][0];
119 simstate.cell_pressure_[cell] = PhaseVec(press);
120 typename Fluid::FluidState state = fluid.computeState(simstate.cell_pressure_[cell], simstate.cell_z_[cell]);
121 fluid_vol_dens = state.total_phase_volume_density_;
122 simstate.cell_z_[cell] *= 1.0/fluid_vol_dens;
123 ++cnt;
124 } while (std::fabs((fluid_vol_dens-1.0)) > 1.0e-8 && cnt < 10);
125
126 }
127 } else {
128 std::cout << "---- Exit - BlackoilSimulator.hpp: No gravity, no fun ... ----" << std::endl;
129 exit(-1);
130 }
131 } else if (param.getDefault("CO2-injection", false)) {
132 CompVec init_water(0.0);
133 // Initially water filled (use Oil-component for water in order
134 // to utilise blackoil mechanisms for brine-co2 interaction)
135 init_water[Fluid::Oil] = 1.0;
136 simstate.cell_z_.resize(grid.numCells());
137 std::fill(simstate.cell_z_.begin(),simstate.cell_z_.end(),init_water);
138
139 double datum_pressure_barsa = param.getDefault<double>("datum_pressure", 200.0);
140 double datum_pressure = Opm::unit::convert::from(datum_pressure_barsa, Opm::unit::barsa);
141 PhaseVec init_p(datum_pressure);
142 simstate.cell_pressure_.resize(grid.numCells(), init_p);
143
144 // Simple initial condition based on "incompressibility"-assumption
145 double zMin = grid.cellCentroid(0)[2];
146 for (int cell = 1; cell < grid.numCells(); ++cell) {
147 if (grid.cellCentroid(cell)[2] < zMin)
148 zMin = grid.cellCentroid(cell)[2];
149 }
150
151 typename Fluid::FluidState state = fluid.computeState(init_p, init_water);
152 simstate.cell_z_[0] *= 1.0/state.total_phase_volume_density_;
153 double density = (init_water*fluid.surfaceDensities())/state.total_phase_volume_density_;
154
155 for (int cell = 0; cell < grid.numCells(); ++cell) {
156 double pressure(datum_pressure + (grid.cellCentroid(cell)[2] - zMin)*gravity[2]*density);
157 simstate.cell_pressure_[cell] = PhaseVec(pressure);
158 state = fluid.computeState(simstate.cell_pressure_[cell], simstate.cell_z_[cell]);
159 simstate.cell_z_[cell] *= 1.0/state.total_phase_volume_density_;
160 }
161 } else {
162 CompVec init_z(0.0);
163 double initial_mixture_gas = param.getDefault("initial_mixture_gas", 0.0);
164 double initial_mixture_oil = param.getDefault("initial_mixture_oil", 1.0);
165 double initial_mixture_water = param.getDefault("initial_mixture_water", 0.0);
166 init_z[Fluid::Water] = initial_mixture_water;
167 init_z[Fluid::Gas] = initial_mixture_gas;
168 init_z[Fluid::Oil] = initial_mixture_oil;
169
170 simstate.cell_z_.resize(grid.numCells(), init_z);
171 OPM_MESSAGE("******* Assuming zero capillary pressures *******");
172 PhaseVec init_p(param.getDefault("initial_pressure", 100.0*Opm::unit::barsa));
173 simstate.cell_pressure_.resize(grid.numCells(), init_p);
174 if (gravity.two_norm() != 0.0) {
175 double ref_gravpot = grid.cellCentroid(0)*gravity;
176 double rho = init_z*fluid.surfaceDensities(); // Assuming incompressible, and constant initial z.
177 for (int cell = 1; cell < grid.numCells(); ++cell) {
178 double press = rho*(grid.cellCentroid(cell)*gravity - ref_gravpot) + simstate.cell_pressure_[0][0];
179 simstate.cell_pressure_[cell] = PhaseVec(press);
180 }
181 }
182 }
183 }
184 };
185
186
187
188
189
190
192 template <class Simulator>
193 class SPE9Initialization : public BlackoilInitialization<Simulator>
194 {
195 public:
196 typedef typename Simulator::State State;
197 typedef typename Simulator::Grid Grid;
198 typedef typename Simulator::Fluid Fluid;
199
200 virtual void init(const Opm::parameter::ParameterGroup& param,
201 const Grid& grid,
202 const Fluid& fluid,
203 typename Grid::Vector gravity,
204 State& simstate)
205 {
206 typedef typename Fluid::CompVec CompVec;
207
208 double zeroDepth = param.getDefault("zero_depth", 2743.2);
209 int nx = param.getDefault<int>("nx", 24);
210 int ny = param.getDefault<int>("ny", 25);
211 int nz = param.getDefault<int>("nz", 15);
212
213 // double datum_depth = param.getDefault<double>("datum_depth", 2753.87) - zeroDepth;
214 double datum_pressure_barsa = param.getDefault<double>("datum_pressure", 248.22);
215 double datum_pressure = Opm::unit::convert::from(datum_pressure_barsa, Opm::unit::barsa);
216 double wo_contact_depth = param.getDefault<double>("wo_contact_depth", 3032.76) - zeroDepth;
217 double go_contact_depth = param.getDefault<double>("go_contact_depth", 2682.24) - zeroDepth;
218
219 double connate_water_saturation = param.getDefault<double>("connate_water_saturation", 0.151090);
220 double residual_oil_saturation = param.getDefault<double>("residual_oil_saturation", 0.118510);
221
222 double initial_mixture_gas = param.getDefault("initial_mixture_gas", 247.43);
223 double initial_mixture_oil = param.getDefault("initial_mixture_oil", 1.0);
224
225 // Initial fluid state
226 CompVec oil_sample(0.0);
227 oil_sample[Fluid::Oil] = initial_mixture_oil;
228 oil_sample[Fluid::Gas] = initial_mixture_gas;
229 CompVec water_sample(0.0);
230 water_sample[Fluid::Water] = 1.0;
231
232 simstate.cell_z_.resize(grid.numCells());
233 simstate.cell_pressure_.resize(grid.numCells());
234
235 // Datum -cell
236 // For now, assume that datum_depth corresponds the centroid of cell 0 (reasonable approx)
237 simstate.cell_pressure_[0] = datum_pressure;
238 typename Fluid::FluidState state = fluid.computeState(simstate.cell_pressure_[0],oil_sample);
239 simstate.cell_z_[0] = oil_sample;
240 simstate.cell_z_[0] *= (1.0-connate_water_saturation)/state.total_phase_volume_density_;
241 state = fluid.computeState(simstate.cell_pressure_[0],water_sample);
242 simstate.cell_z_[0][Fluid::Water] = water_sample[Fluid::Water];
243 simstate.cell_z_[0][Fluid::Water] *= connate_water_saturation/state.total_phase_volume_density_;
244 // Rest of the cells -- NOTE: Assume uniform cell properties in y-direction
245 for (int i=0; i<nx; ++i) {
246 int k0=i*nz;
247 for (int k=0; k<nz; ++k) {
248 int kk=k0+k;
249 if (i>0 && k==0) {
250 computeCellState(grid, fluid, gravity,
251 kk, kk-nz, wo_contact_depth, go_contact_depth, connate_water_saturation,
252 residual_oil_saturation, simstate);
253 } else if (k>0) {
254 computeCellState(grid, fluid, gravity,
255 kk, kk-1, wo_contact_depth, go_contact_depth, connate_water_saturation,
256 residual_oil_saturation, simstate);
257 }
258 // Copy cell properties to y-layers
259 for (int j=1; j<ny; ++j) {
260 int jj = j*nx*nz + kk;
261 simstate.cell_z_[jj] = simstate.cell_z_[kk];
262 simstate.cell_pressure_[jj] = simstate.cell_pressure_[kk];
263 }
264 }
265 }
266 }
267
268
269
270 bool computeCellState(const Grid& grid,
271 const Fluid& fluid,
272 typename Grid::Vector gravity,
273 int iCell,
274 int iRef,
275 double wo_contact_depth,
276 double go_contact_depth,
277 double connate_water_saturation,
278 double residual_oil_saturation,
279 State& simstate)
280 {
281 typedef typename Fluid::PhaseVec PhaseVec;
282
283 const int maxCnt = 30;
284 const double eps = 1.0e-8;
285 simstate.cell_z_[iCell] = simstate.cell_z_[iRef];
286
287 bool below_wo_contact = false;
288 if (grid.cellCentroid(iCell)[2] > wo_contact_depth)
289 below_wo_contact = true;
290
291 double gZ = (grid.cellCentroid(iCell) - grid.cellCentroid(iRef))*gravity;
292 double fluid_vol_dens;
293 int cnt =0;
294 do {
295 double rho = 0.5*(simstate.cell_z_[iCell]*fluid.surfaceDensities()
296 + simstate.cell_z_[iRef]*fluid.surfaceDensities());
297 double press = rho*gZ + simstate.cell_pressure_[iRef][0];
298 simstate.cell_pressure_[iCell] = PhaseVec(press);
299 typename Fluid::FluidState state = fluid.computeState(simstate.cell_pressure_[iCell], simstate.cell_z_[iCell]);
300 fluid_vol_dens = state.total_phase_volume_density_;
301 double oil_vol_dens = state.phase_volume_density_[Fluid::Liquid]
302 + state.phase_volume_density_[Fluid::Vapour];
303 double wat_vol_dens = state.phase_volume_density_[Fluid::Aqua];
304 if (below_wo_contact) {
305 simstate.cell_z_[iCell][Fluid::Oil] *= residual_oil_saturation/oil_vol_dens;
306 simstate.cell_z_[iCell][Fluid::Gas] *= residual_oil_saturation/oil_vol_dens;
307 simstate.cell_z_[iCell][Fluid::Water] *= (1.0-residual_oil_saturation)/wat_vol_dens;
308 } else {
309 simstate.cell_z_[iCell][Fluid::Oil] *= (1.0-connate_water_saturation)/oil_vol_dens;
310 simstate.cell_z_[iCell][Fluid::Gas] *= (1.0-connate_water_saturation)/oil_vol_dens;
311 simstate.cell_z_[iCell][Fluid::Water] *= connate_water_saturation/wat_vol_dens;
312 }
313 ++cnt;
314 } while (std::fabs(fluid_vol_dens-1.0) > eps && cnt < maxCnt);
315
316 if (cnt == maxCnt) {
317 std::cout << "z_cell_[" << iCell << "]: " << simstate.cell_z_[iCell]
318 << " pressure: " << simstate.cell_pressure_[iCell][Fluid::Liquid]
319 << " cnt: " << cnt
320 << " eps: " << std::fabs(fluid_vol_dens-1.0) << std::endl;
321 }
322
323 return (cnt < maxCnt);
324 }
325
326 };
327
328} // namespace Opm
329
330#endif // OPM_BLACKOILINITIALIZATION_HEADER_INCLUDED
Opm::BasicInitialization
Initialize basic cases.
Definition: BlackoilInitialization.hpp:53
Opm::BasicInitialization::State
Simulator::State State
Definition: BlackoilInitialization.hpp:55
Opm::BasicInitialization::init
virtual void init(const Opm::parameter::ParameterGroup &param, const Grid &grid, const Fluid &fluid, typename Grid::Vector gravity, State &simstate)
Definition: BlackoilInitialization.hpp:59
Opm::BasicInitialization::Fluid
Simulator::Fluid Fluid
Definition: BlackoilInitialization.hpp:57
Opm::BasicInitialization::Grid
Simulator::Grid Grid
Definition: BlackoilInitialization.hpp:56
Opm::BlackoilInitialization
Base class for initialization codes.
Definition: BlackoilInitialization.hpp:35
Opm::BlackoilInitialization::State
Simulator::State State
Definition: BlackoilInitialization.hpp:37
Opm::BlackoilInitialization::init
virtual void init(const Opm::parameter::ParameterGroup &param, const Grid &grid, const Fluid &fluid, typename Grid::Vector gravity, State &simstate)=0
Opm::BlackoilInitialization::Fluid
Simulator::Fluid Fluid
Definition: BlackoilInitialization.hpp:39
Opm::BlackoilInitialization::Grid
Simulator::Grid Grid
Definition: BlackoilInitialization.hpp:38
Opm::SPE9Initialization
Initialize SPE9 type case.
Definition: BlackoilInitialization.hpp:194
Opm::SPE9Initialization::init
virtual void init(const Opm::parameter::ParameterGroup &param, const Grid &grid, const Fluid &fluid, typename Grid::Vector gravity, State &simstate)
Definition: BlackoilInitialization.hpp:200
Opm::SPE9Initialization::State
Simulator::State State
Definition: BlackoilInitialization.hpp:196
Opm::SPE9Initialization::computeCellState
bool computeCellState(const Grid &grid, const Fluid &fluid, typename Grid::Vector gravity, int iCell, int iRef, double wo_contact_depth, double go_contact_depth, double connate_water_saturation, double residual_oil_saturation, State &simstate)
Definition: BlackoilInitialization.hpp:270
Opm::SPE9Initialization::Grid
Simulator::Grid Grid
Definition: BlackoilInitialization.hpp:197
Opm::SPE9Initialization::Fluid
Simulator::Fluid Fluid
Definition: BlackoilInitialization.hpp:198
Opm
Definition: BlackoilFluid.hpp:32