57 enum { dimWorld = GridView::dimensionworld };
58 enum { numPhases = FluidSystem::numPhases };
59 static constexpr int numEq = BlackoilIndices::numEq;
61 using Eval = DenseAd::Evaluation<Scalar, numEq>;
62 using Toolbox = MathToolbox<Eval>;
64 using typename AquiferInterface<TypeTag>::RateVector;
68 const Simulator& simulator)
71 , cumulative_flux_(0.0)
72 , init_pressure_ (aquifer.numCells(), 0.0)
74 this->cell_to_aquifer_cell_idx_.resize(this->simulator_.gridView().size(0), -1);
76 auto aquifer_on_process =
false;
77 for (std::size_t idx = 0; idx < aquifer.numCells(); ++idx) {
78 const auto* cell = aquifer.getCellPrt(idx);
81 const int compressed_idx = simulator.vanguard().compressedIndexForInterior(cell->global_index);
82 if (compressed_idx >= 0) {
83 this->cell_to_aquifer_cell_idx_[compressed_idx] = idx;
84 aquifer_on_process =
true;
88 if (aquifer_on_process) {
89 this->checkConnectsToReservoir();
96 result.flux_rate_ = 1.0;
97 result.cumulative_flux_ = 2.0;
98 result.init_pressure_ = {3.0, 4.0};
99 result.pressure_ = 5.0;
104 void initFromRestart(
const data::Aquifers& aquiferSoln)
override
106 auto xaqPos = aquiferSoln.find(this->aquiferID());
107 if (xaqPos == aquiferSoln.end())
110 if (this->connects_to_reservoir_) {
111 this->cumulative_flux_ = xaqPos->second.volume;
114 if (
const auto* aqData = xaqPos->second.typeData.template get<data::AquiferType::Numerical>();
117 this->init_pressure_.resize(aqData->initPressure.size());
118 std::copy(aqData->initPressure.begin(),
119 aqData->initPressure.end(),
120 this->init_pressure_.begin());
123 this->solution_set_from_restart_ =
true;
126 void beginTimeStep()
override {}
127 void addToSource(RateVector&,
const unsigned,
const unsigned)
override {}
129 void endTimeStep()
override
131 this->pressure_ = this->calculateAquiferPressure();
132 this->flux_rate_ = this->calculateAquiferFluxRate();
133 this->cumulative_flux_ += this->flux_rate_ * this->simulator_.timeStepSize();
136 data::AquiferData aquiferData()
const override
138 data::AquiferData data;
139 data.aquiferID = this->aquiferID();
140 data.pressure = this->pressure_;
141 data.fluxRate = this->flux_rate_;
142 data.volume = this->cumulative_flux_;
144 auto* aquNum = data.typeData.template create<data::AquiferType::Numerical>();
145 aquNum->initPressure.resize(this->init_pressure_.size());
146 std::copy(this->init_pressure_.begin(),
147 this->init_pressure_.end(),
148 aquNum->initPressure.begin());
153 void initialSolutionApplied()
override
155 if (this->solution_set_from_restart_) {
159 this->pressure_ = this->calculateAquiferPressure(this->init_pressure_);
160 this->flux_rate_ = 0.;
161 this->cumulative_flux_ = 0.;
164 void computeFaceAreaFraction(
const std::vector<Scalar>& )
override
167 Scalar totalFaceArea()
const override
172 template<
class Serializer>
173 void serializeOp(Serializer& serializer)
175 serializer(flux_rate_);
176 serializer(cumulative_flux_);
177 serializer(init_pressure_);
178 serializer(pressure_);
183 return this->flux_rate_ == rhs.flux_rate_ &&
184 this->cumulative_flux_ == rhs.cumulative_flux_ &&
185 this->init_pressure_ == rhs.init_pressure_ &&
186 this->pressure_ == rhs.pressure_;
189 Scalar cumulativeFlux()
const
191 return this->cumulative_flux_;
195 void checkConnectsToReservoir()
197 ElementContext elem_ctx(this->simulator_);
198 auto elemIt = std::find_if(this->simulator_.gridView().template begin</*codim=*/0>(),
199 this->simulator_.gridView().template end</*codim=*/0>(),
200 [&elem_ctx,
this](
const auto& elem) ->
bool
202 elem_ctx.updateStencil(elem);
204 const auto cell_index = elem_ctx
205 .globalSpaceIndex(0, 0);
207 return this->cell_to_aquifer_cell_idx_[cell_index] == 0;
210 assert ((elemIt != this->simulator_.gridView().template end</*codim=*/0>())
211 &&
"Internal error locating numerical aquifer's connecting cell");
213 this->connects_to_reservoir_ =
214 elemIt->partitionType() == Dune::InteriorEntity;
217 Scalar calculateAquiferPressure()
const
219 auto capture = std::vector<Scalar>(this->init_pressure_.size(), 0.0);
220 return this->calculateAquiferPressure(capture);
223 Scalar calculateAquiferPressure(std::vector<Scalar>& cell_pressure)
const
225 Scalar sum_pressure_watervolume = 0.;
226 Scalar sum_watervolume = 0.;
228 ElementContext elem_ctx(this->simulator_);
229 const auto& gridView = this->simulator_.gridView();
230 OPM_BEGIN_PARALLEL_TRY_CATCH();
232 for (
const auto& elem : elements(gridView, Dune::Partitions::interior)) {
233 elem_ctx.updatePrimaryStencil(elem);
235 const std::size_t cell_index = elem_ctx.globalSpaceIndex(0, 0);
236 const int idx = this->cell_to_aquifer_cell_idx_[cell_index];
241 elem_ctx.updatePrimaryIntensiveQuantities(0);
242 const auto& iq0 = elem_ctx.intensiveQuantities(0, 0);
243 const auto& fs = iq0.fluidState();
248 const Scalar water_saturation = fs.saturation(this->phaseIdx_()).value();
249 const Scalar porosity = iq0.porosity().value();
250 const Scalar volume = elem_ctx.dofTotalVolume(0, 0);
252 const Scalar water_pressure_reservoir = fs.pressure(this->phaseIdx_()).value();
253 const Scalar water_volume = volume * porosity * water_saturation;
254 sum_pressure_watervolume += water_volume * water_pressure_reservoir;
255 sum_watervolume += water_volume;
257 cell_pressure[idx] = water_pressure_reservoir;
259 OPM_END_PARALLEL_TRY_CATCH(
"AquiferNumerical::calculateAquiferPressure() failed: ",
260 this->simulator_.vanguard().grid().comm());
261 const auto& comm = this->simulator_.vanguard().grid().comm();
262 comm.sum(&sum_pressure_watervolume, 1);
263 comm.sum(&sum_watervolume, 1);
266 comm.sum(cell_pressure.data(), cell_pressure.size());
268 return sum_pressure_watervolume / sum_watervolume;
271 template <
class ElemCtx>
272 Scalar getWaterFlux(
const ElemCtx& elem_ctx,
unsigned face_idx)
const
274 const auto& exQuants = elem_ctx.extensiveQuantities(face_idx, 0);
275 const Scalar water_flux = Toolbox::value(exQuants.volumeFlux(this->phaseIdx_()));
279 Scalar calculateAquiferFluxRate()
const
281 Scalar aquifer_flux = 0.0;
283 if (! this->connects_to_reservoir_) {
287 ElementContext elem_ctx(this->simulator_);
288 const auto& gridView = this->simulator_.gridView();
289 for (
const auto& elem : elements(gridView, Dune::Partitions::interior)) {
290 elem_ctx.updatePrimaryStencil(elem);
291 const std::size_t cell_index = elem_ctx.globalSpaceIndex(0, 0);
292 const int idx = this->cell_to_aquifer_cell_idx_[cell_index];
298 elem_ctx.updateStencil(elem);
299 const std::size_t num_interior_faces = elem_ctx.numInteriorFaces( 0);
300 const auto& stencil = elem_ctx.stencil(0);
301 elem_ctx.updateAllIntensiveQuantities();
302 elem_ctx.updateAllExtensiveQuantities();
304 for (std::size_t face_idx = 0; face_idx < num_interior_faces; ++face_idx) {
305 const auto& face = stencil.interiorFace(face_idx);
307 const std::size_t i = face.interiorIndex();
308 const std::size_t j = face.exteriorIndex();
311 const std::size_t J = stencil.globalSpaceIndex(j);
313 assert(stencil.globalSpaceIndex(i) == cell_index);
317 if (this->cell_to_aquifer_cell_idx_[J] > 0) {
321 const Scalar water_flux = getWaterFlux(elem_ctx,face_idx);
322 const std::size_t up_id = water_flux >= 0.0 ? i : j;
323 const auto& intQuantsIn = elem_ctx.intensiveQuantities(up_id, 0);
324 const Scalar invB = Toolbox::value(intQuantsIn.fluidState().invB(this->phaseIdx_()));
325 const Scalar face_area = face.area();
326 aquifer_flux += water_flux * invB * face_area;
337 Scalar cumulative_flux_;
338 std::vector<Scalar> init_pressure_{};
340 bool solution_set_from_restart_ {
false};
341 bool connects_to_reservoir_ {
false};
344 std::vector<int> cell_to_aquifer_cell_idx_;