diff --git a/opm/geomech/Fracture.cpp b/opm/geomech/Fracture.cpp index 24b665b..c49094f 100644 --- a/opm/geomech/Fracture.cpp +++ b/opm/geomech/Fracture.cpp @@ -1,30 +1,39 @@ #include "config.h" -#include "opm/geomech/GridStretcher.hpp" -#include + #include -#include -#include -#include -#include + #include // needed for printSparseMatrix?? +#include +#include #include // needed for printSparseMatrix?? + +#include + +#include +#include + +#include +#include +#include + +#include #include #include #include #include -#include -#include +#include +#include -#include +#include namespace { /** * @brief Computes the target expansion for a fracture. * - * This function calculates the target expansion based on the given parameters: - * the target stress intensity factor (K1_target), the fracture aperture, - * Young's modulus (E), and Poisson's ratio (nu). + * This function calculates the target expansion based on the given + * parameters: the target stress intensity factor (K1_target), the fracture + * aperture, Young's modulus (E), and Poisson's ratio (nu). * * @param K1_target The target stress intensity factor. * @param aperture The fracture aperture. @@ -33,55 +42,65 @@ namespace { * @return The computed target expansion. */ double compute_target_expansion(const double K1_target, - const double aperture, - const double E, // young - const double nu) // poisson + const double aperture, + const double E, // young + const double nu) // poisson { - const double mu = E / (2 * (1+nu)); // shear modulus - const double fac = mu * std::sqrt(M_PI) / - (2 * (1-nu) * 1.834); - return pow(fac * aperture / K1_target, 2); -}; - -}; // end anonymous namespace + const double mu = E / (2 * (1 + nu)); // shear modulus + const double fac = mu * std::sqrt(M_PI) / + (2 * (1 - nu) * 1.834); + + return std::pow(fac * aperture / K1_target, 2); +} + +} // end anonymous namespace namespace Opm { -void -Fracture::init(std::string well, - int perf, - int well_cell, - Fracture::Point3D origo, - Fracture::Point3D normal, - Opm::PropertyTree prm - ) +void Fracture::init(const std::string well, + const int perf, + const int well_cell, + const Fracture::Point3D origo, + const Fracture::Point3D normal, + Opm::PropertyTree prm) { prm_ = prm; wellinfo_ = WellInfo({well, perf, well_cell}); origo_ = origo; + axis_[2] = normal; - axis_[0] = Point3D({std::copysign(normal[2], normal[0]), - std::copysign(normal[2], normal[1]), - -std::copysign(std::abs(normal[0]) + std::abs(normal[1]), normal[2])}); + axis_[0] = Point3D { + std::copysign(normal[2], normal[0]), + std::copysign(normal[2], normal[1]), + -std::copysign(std::abs(normal[0]) + std::abs(normal[1]), normal[2]) + }; + axis_[1] = crossProduct(axis_[2], axis_[0]); - double init_scale = prm_.get("config.axis_scale"); + + const double init_scale = prm_.get("config.axis_scale"); for (int i = 0; i < 3; ++i) { axis_[i] /= axis_[i].two_norm(); axis_[i] *= init_scale; } + layers_ = 0; nlinear_ = 0; + initFracture(); - grow(1, 0); + this->grow(1, 0); + nlinear_ = layers_; + //grow(4, 1); - grid_->grow(); - grid_->postGrow(); + this->grid_->grow(); + this->grid_->postGrow(); - size_t nc = grid_->leafGridView().size(0); + const size_t nc = grid_->leafGridView().size(0); reservoir_cells_ = std::vector(nc, -1); - fracture_pressure_.resize(nc); fracture_pressure_ = 1e5; + fracture_pressure_.resize(nc); + fracture_pressure_ = 1e5; + this->initFractureWidth(); //fracture_width_.resize(nc); fracture_width_ = 1e-3; // @@ -91,77 +110,76 @@ Fracture::init(std::string well, //setGridStretcher(); this->resetWriters(); - - } - void Fracture::setGridStretcher() { grid_stretcher_ = std::unique_ptr(new GridStretcher(*grid_)); } - -void Fracture::resetWriters(){ +void Fracture::resetWriters() +{ // nead to be reseat if grid is changed ?? vtkwriter_ = - std::make_unique>(grid_->leafGridView(), - Dune::VTK::nonconforming); + std::make_unique> + (grid_->leafGridView(), Dune::VTK::nonconforming); + + const std::string outputdir = prm_.get("outputdir"); + const std::string simName = prm_.get("casename") + this->name(); + const std::string multiFileName = ""; - std::string outputdir = prm_.get("outputdir"); - std::string simName = prm_.get("casename") + this->name(); - std::string multiFileName = ""; vtkmultiwriter_ = - std::make_unique< Opm::VtkMultiWriter >( - /*async*/ false, - grid_->leafGridView(), - outputdir, - simName, - multiFileName ); + std::make_unique> + (/*async*/ false, + grid_->leafGridView(), + outputdir, + simName, + multiFileName); } -void Fracture::setupPressureSolver(){ +void Fracture::setupPressureSolver() +{ Opm::FlowLinearSolverParameters p; p.linsolver_ = prm_.get("pressuresolver"); prmpressure_ = Opm::setupPropertyTree(p, true, true); - { - std::size_t pressureIndex; // Dummy - const std::function weightsCalculator; // Dummy - auto pressure_operator = std::make_unique(*pressure_matrix_); - // using FlexibleSolverType = Dune::FlexibleSolver; - pressure_operator_ = std::move(pressure_operator); - auto psolver - = std::make_unique(*pressure_operator_, prmpressure_, weightsCalculator, pressureIndex); - - pressure_solver_ = std::move(psolver); - } + + this->pressure_operator_ = + std::make_unique(*pressure_matrix_); + + const std::size_t pressureIndex{}; // Dummy + const std::function weightsCalculator{}; // Dummy + this->pressure_solver_ = std::make_unique + (*pressure_operator_, prmpressure_, weightsCalculator, pressureIndex); } + /** * @brief Removes cells from the grid if they are out side reservoir. * * This function performs the following steps: * 1. Copies the current fracture width data to a persistent container. - * 2. Iterates over all elements in the grid's leaf view and removes elements where the reservoir cell index is negative. + * 2. Iterates over all elements in the grid's leaf view and removes + * elements where the reservoir cell index is negative. * 3. Grows the grid and performs post-growth operations. * 4. Resizes the fracture width array to match the new grid size. - * 5. Copies the fracture width data back from the persistent container to the resized array. + * 5. Copies the fracture width data back from the persistent container to + * the resized array. * 6. Resets the writers associated with the Fracture object. */ -void Fracture::removeCells(){ +void Fracture::removeCells() +{ // copy all to presistent container ElementMapper mapper(grid_->leafGridView(), Dune::mcmgElementLayout()); // used id sets interally - Dune::PersistentContainer fracture_width(*grid_,0); + + Dune::PersistentContainer fracture_width(*grid_, 0); fracture_width.resize(); - for(const auto& elem: elements(grid_->leafGridView())){ - size_t eIdx = mapper.index(elem); - fracture_width[elem] = fracture_width_[eIdx]; + for (const auto& elem : elements(grid_->leafGridView())) { + fracture_width[elem] = fracture_width_[mapper.index(elem)]; } //const auto& indexSet = foamGridLeafView.indexSet();// for indices //const auto& indexSet = grid.localIdSet();// presitent numbering - for(const auto& elem: elements(grid_->leafGridView())){ - size_t eIdx = mapper.index(elem); - if(reservoir_cells_[eIdx] < 0){ + for (const auto& elem : elements(grid_->leafGridView())) { + if (reservoir_cells_[mapper.index(elem)] < 0) { grid_->removeElement(elem); } } @@ -170,60 +188,56 @@ void Fracture::removeCells(){ grid_->postGrow(); // resize array - fracture_width_.resize(grid_->leafGridView().size(0)); // copy back from presistent contatiner - for(const auto& elem: elements(grid_->leafGridView())){ - size_t eIdx = mapper.index(elem); - fracture_width_[eIdx] = fracture_width[elem]; + fracture_width_.resize(grid_->leafGridView().size(0)); + for (const auto& elem : elements(grid_->leafGridView())) { + fracture_width_[mapper.index(elem)] = fracture_width[elem]; } + this->resetWriters(); } -Dune::BlockVector> Fracture::all_slips() const{ +Dune::BlockVector> Fracture::all_slips() const +{ Dune::BlockVector> slips(grid_->leafGridView().size(0)); - slips = 0; + slips = 0.0; + ElementMapper mapper(grid_->leafGridView(), Dune::mcmgElementLayout()); - for(const auto& elem: Dune::elements(grid_->leafGridView())){ - size_t eIdx = mapper.index(elem); + for (const auto& elem : Dune::elements(grid_->leafGridView())) { + const auto eIdx = mapper.index(elem); //only normal slip for now slips[eIdx][0] = fracture_width_[eIdx]; } + return slips; } -Dune::FieldVector Fracture::disp(Dune::FieldVector obs) const{ - auto slips = all_slips(); - Dune::FieldVector disp = ddm::disp(obs, slips, *grid_, E_, nu_); - return disp; +Dune::FieldVector Fracture::disp(const Dune::FieldVector& obs) const +{ + return ddm::disp(obs, this->all_slips(), *grid_, E_, nu_); } -Dune::FieldVector Fracture::strain(Dune::FieldVector obs) const{ +Dune::FieldVector Fracture::strain(const Dune::FieldVector& obs) const +{ // for now use full slip in interface even we only calculate normal slip - Dune::BlockVector> slips = all_slips(); - Dune::FieldVector strain = ddm::strain(obs, slips, *grid_, E_, nu_); - return strain; + return ddm::strain(obs, this->all_slips(), *grid_, E_, nu_); } -Dune::FieldVector Fracture::stress(Dune::FieldVector obs) const{ - // for now use full slip in interface even we only calculate normal slip - Dune::BlockVector> slips = all_slips(); - Dune::FieldVector strain = ddm::strain(obs, slips, *grid_, E_, nu_); - Dune::FieldVector stress = ddm::strainToStress(E_, nu_, strain); - return stress; +Dune::FieldVector Fracture::stress(const Dune::FieldVector& obs) const +{ + return ddm::strainToStress(E_, nu_, this->strain(obs)); } -std::string -Fracture::name() const +std::string Fracture::name() const { - std::string name = "Fracure_on_" + wellinfo_.name + "_perf_" + std::to_string(wellinfo_.perf) + "_nr" ; - return name; + return "Fracure_on_" + wellinfo_.name + "_perf_" + std::to_string(wellinfo_.perf) + "_nr" ; } -void -Fracture::initFracture() + +void Fracture::initFracture() { Dune::GridFactory factory; // Dune::FoamGrid<2,3>> factory; - size_t N = 6; - double radius = 1; + const std::size_t N = 6; + const double radius = 1.0; std::vector inner_indices; std::vector> vertices; { @@ -232,58 +246,65 @@ Fracture::initFracture() vertices.push_back(vertex); } - for (size_t i = 0; i < N; ++i) { - + for (std::size_t i = 0; i < N; ++i) { inner_indices.push_back(i + 1); - double theta = (i * 2 * M_PI) / N; - double x = radius * cos(theta); - double y = radius * sin(theta); - Point3D vertex = surfaceMap(x, y); - vertices.push_back(vertex); + + const double theta = (i * 2 * M_PI) / N; + const double x = radius * cos(theta); + const double y = radius * sin(theta); + + vertices.push_back(surfaceMap(x, y)); // assume the first 6 cells has source well_source_.push_back(i); } - for (size_t i = 0; i < vertices.size(); i++) { - factory.insertVertex(vertices[i]); + + for (const auto& vertex : vertices) { + factory.insertVertex(vertex); } + std::vector> cornerIDs; - for (size_t i = 0; i < N; ++i) { + for (std::size_t i = 0; i < N; ++i) { unsigned int next = inner_indices[(i + 1) % N]; std::vector cornerID({unsigned(0), unsigned(inner_indices[i]), next}); cornerIDs.push_back(cornerID); } - for (size_t i = 0; i < N; ++i) { + + for (std::size_t i = 0; i < N; ++i) { factory.insertElement(Dune::GeometryTypes::simplex(2), cornerIDs[i]); // std::move(mapping)); } + out_indices_ = inner_indices; grid_ = factory.createGrid(); // grid_ = factory.createGrid(); } -std::vector Fracture::stressIntensityK1() const{ - size_t nc = grid_->leafGridView().size(0); - std::vector stressIntensityK1(nc, std::nan("0")); - ElementMapper mapper(grid_->leafGridView(), Dune::mcmgElementLayout()); - for(const auto& elem: elements(grid_->leafGridView())){ - //bool isboundary = false; - for (auto& is : Dune::intersections(grid_->leafGridView(),elem)) { - if (is.boundary()) { - int nIdx = mapper.index(elem); - auto isCenter = is.geometry().center(); - auto elCenter = elem.geometry().center(); - auto vecC = isCenter-elCenter; - auto distC = vecC.two_norm(); - double K1 = ddm::fractureK1(distC, fracture_width_[nIdx], E_, nu_); - stressIntensityK1[nIdx] = K1; - } +std::vector Fracture::stressIntensityK1() const +{ + auto stressIntensityK1 = std::vector + (grid_->leafGridView().size(0), std::nan("0")); + + ElementMapper mapper(grid_->leafGridView(), Dune::mcmgElementLayout()); + for (const auto& elem : elements(grid_->leafGridView())) { + const auto nIdx = mapper.index(elem); + + const auto elCenter = elem.geometry().center(); + const auto fracWidth = this->fracture_width_[nIdx]; + + auto& K1 = stressIntensityK1[nIdx]; + for (auto& is : Dune::intersections(grid_->leafGridView(), elem)) { + if (! is.boundary()) { + continue; } + const auto vecC = is.geometry().center() - elCenter; + K1 = ddm::fractureK1(vecC.two_norm(), fracWidth, E_, nu_); } - return stressIntensityK1; } -void -Fracture::write(int reportStep) const + return stressIntensityK1; +} + +void Fracture::write(int reportStep) const { std::vector K1;//need to be in scope until written if (reservoir_cells_.size() > 0) { @@ -362,7 +383,7 @@ void Fracture::writemulti(double time) const // maybe slow well_index[i] = wellIndMap[reservoir_cells_[i]]; } - + vtkmultiwriter_->beginWrite(time); @@ -405,8 +426,7 @@ void Fracture::writemulti(double time) const vtkmultiwriter_->endWrite(); }; -void -Fracture::grow(int layers, int method) +void Fracture::grow(int layers, int method) { while (layers_ < layers) { std::vector inner_indices = out_indices_; @@ -421,91 +441,124 @@ Fracture::grow(int layers, int method) inner_indices = out_indices_; } } -void -Fracture::insertLinear(const std::vector& inner_indices) + +void Fracture::insertLinear(const std::vector& inner_indices) { - double radius = 1.0; - size_t N = inner_indices.size(); - for (size_t i = 0; i < N; ++i) { + const double radius = 1.0; + const auto N = inner_indices.size(); + + for (auto i = 0*N; i < N; ++i) { double theta = (i * 2 * M_PI) / (N); theta += (layers_ - nlinear_) * 0.5 * (2 * M_PI / N); - double out_radius = radius + layers_ + 1; - double x = out_radius * cos(theta); - double y = out_radius * sin(theta); - Point3D vertex = surfaceMap(x, y); - int new_ind = grid_->insertVertex(vertex); + + const double out_radius = radius + layers_ + 1; + const double x = out_radius * std::cos(theta); + const double y = out_radius * std::sin(theta); + + const Point3D vertex = surfaceMap(x, y); + const int new_ind = grid_->insertVertex(vertex); + out_indices_.push_back(new_ind); } - for (size_t i = 0; i < N; ++i) { + + for (auto i = 0*N; i < N; ++i) { { - std::vector cornerID( - {inner_indices[i % N], out_indices_[(i) % (N)], inner_indices[(i + 1) % (N)]}); + const std::vector cornerID { + inner_indices[i % N], + out_indices_[(i) % (N)], + inner_indices[(i + 1) % (N)] + }; + grid_->insertElement(Dune::GeometryTypes::simplex(2), cornerID); } { - std::vector cornerID( - {inner_indices[(i + 1) % N], out_indices_[(i) % (N)], out_indices_[(i + 1) % (N)]}); + const std::vector cornerID { + inner_indices[(i + 1) % N], + out_indices_[(i) % (N)], + out_indices_[(i + 1) % (N)] + }; + grid_->insertElement(Dune::GeometryTypes::simplex(2), cornerID); } } } -void -Fracture::insertExp(const std::vector& inner_indices) +void Fracture::insertExp(const std::vector& inner_indices) { - size_t N = inner_indices.size(); - double radius = 1.0; - - for (size_t i = 0; i < N * 2; ++i) { - double theta = (i * 2 * M_PI) / (N * 2); - double out_radius = radius + layers_ + 1; - double x = out_radius * cos(theta); - double y = out_radius * sin(theta); - Point3D vertex = surfaceMap(x, y); - int new_ind = grid_->insertVertex(vertex); + const auto N = inner_indices.size(); + const double radius = 1.0; + + for (auto i = 0*N; i < N * 2; ++i) { + const double theta = (i * 2 * M_PI) / (N * 2); + + const double out_radius = radius + layers_ + 1; + const double x = out_radius * std::cos(theta); + const double y = out_radius * std::sin(theta); + const Point3D vertex = surfaceMap(x, y); + + const int new_ind = grid_->insertVertex(vertex); out_indices_.push_back(new_ind); } - - for (size_t i = 0; i < N; ++i) { + for (auto i = 0*N; i < N; ++i) { { - std::vector cornerID( - {inner_indices[i], out_indices_[(2 * i) % (N * 2)], out_indices_[(2 * i + 1) % (N * 2)]}); + const std::vector cornerID { + inner_indices[i], + out_indices_[(2 * i) % (N * 2)], + out_indices_[(2 * i + 1) % (N * 2)] + }; + grid_->insertElement(Dune::GeometryTypes::simplex(2), cornerID); } + { - std::vector cornerID( - {inner_indices[i], out_indices_[(2 * i + 1) % (N * 2)], inner_indices[(i + 1) % N]}); + const std::vector cornerID { + inner_indices[i], + out_indices_[(2 * i + 1) % (N * 2)], + inner_indices[(i + 1) % N] + }; + grid_->insertElement(Dune::GeometryTypes::simplex(2), cornerID); } + { - std::vector cornerID( - {inner_indices[(i + 1) % N], out_indices_[(2 * i + 1) % (N * 2)], out_indices_[(2 * i + 2) % (N * 2)]}); + const std::vector cornerID { + inner_indices[(i + 1) % N], + out_indices_[(2 * i + 1) % (N * 2)], + out_indices_[(2 * i + 2) % (N * 2)] + }; + grid_->insertElement(Dune::GeometryTypes::simplex(2), cornerID); } } } -Dune::FieldVector -Fracture::surfaceMap(double x, double y) +Dune::FieldVector Fracture::surfaceMap(double x, double y) { Point3D vec(0.0); + vec += x * axis_[0]; vec += y * axis_[1]; + // for (int i = 0; i < 2; ++i) { // } vec += origo_; + return vec; } + template void -Fracture::updateReservoirCells(const external::cvf::ref& cellSearchTree,const Grid3D& grid3D) +Fracture::updateReservoirCells(const external::cvf::ref& cellSearchTree, + const Grid3D& grid3D) { cell_normals_.resize(grid_->leafGridView().size(0)); reservoir_cells_.resize(grid_->leafGridView().size(0)); + using GridView = typename Grid::LeafGridView; using ElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper; + ElementMapper elemMapper(grid_->leafGridView(), Dune::mcmgElementLayout()); int tri_divide = 0; int tri_outside = 0; @@ -562,13 +615,13 @@ Fracture::updateReservoirCells(const external::cvf::refname() << " : " << tri_divide << " triangles should be devided" - << std::endl; - std::cout << "For Fracture : " << this->name() << " : " << tri_outside << " triangles outside" << std::endl; - std::cout << "Total triangles: " << grid_->leafGridView().size(0) << std::endl; - auto it = std::find(reservoir_cells_.begin(),reservoir_cells_.end(),-1); - auto extended_fractures = prm_.get("extended_fractures"); - if(!(it == reservoir_cells_.end()) && !(extended_fractures) ){ + + std::cout << "For Fracture : " << this->name() << " : " << tri_divide << " triangles should be divided\n" + << "For Fracture : " << this->name() << " : " << tri_outside << " triangles outside\n" + << "Total triangles: " << grid_->leafGridView().size(0) << std::endl; + + auto it = std::find(reservoir_cells_.begin(), reservoir_cells_.end(), -1); + if ((it != reservoir_cells_.end()) && !prm_.get("extended_fractures")) { std::cout << "Remove fracture outside of model" << std::endl; // remove fracture outside of model this->removeCells(); @@ -576,17 +629,16 @@ Fracture::updateReservoirCells(const external::cvf::ref& gptr) { grid_ = std::move(gptr); layers_ = 0; nlinear_ = 0; - size_t nc = grid_->leafGridView().size(0); + const size_t nc = grid_->leafGridView().size(0); fracture_pressure_.resize(nc); fracture_pressure_ = 1e5; reservoir_cells_ = std::vector(nc, -1); - + updateReservoirProperties(); updateGridDiscretizations(); // assemble mechanics and pressure matrices setGridStretcher(); @@ -594,8 +646,7 @@ void Fracture::setFractureGrid(std::unique_ptr& gptr) this->resetWriters(); } -void -Fracture::updateReservoirProperties() +void Fracture::updateReservoirProperties() { // updater for standalone test double perm = prm_.get("reservoir.perm"); @@ -614,10 +665,10 @@ Fracture::updateReservoirProperties() E_ = 1e9; this->initFractureWidth(); } -void -Fracture::solve() + +void Fracture::solve() { - std::cout << "Solve Fracture Pressure" << std::endl; + std::cout << "Solve Fracture Pressure" << std::endl; std::string method = prm_.get("solver.method"); if(method == "nothing"){ }else if(method == "simple"){ @@ -658,7 +709,7 @@ Fracture::solve() fracture_pressure_ = 0.0; // the following values are only relevant if propagation is requested - const int max_iter = 100; + const int max_iter = 100; const double diameter = 2; // @@ compute this from boundary nodes const double tol = 1e-3 * diameter; //1e-5; //1e-5; // @@ const double efac = 2; // @@ heuristic @@ -672,7 +723,7 @@ Fracture::solve() std::vector displacements(N, {0, 0, 0}); std::ofstream k1log("k1log"); // @@ - std::ofstream texplog("texplog"); //@@ + std::ofstream texplog("texplog"); //@@ std::ofstream distlog("distlog"); //@@ int count = 0; while (true) { @@ -705,8 +756,8 @@ Fracture::solve() K1.push_back(K1_not_nan[i]); // loop over cells, determine how much they should be expanded or contracted - for (size_t i = 0; i != N; ++i) - cell_disp[i] = + for (size_t i = 0; i != N; ++i) + cell_disp[i] = efac * (compute_target_expansion(K1max, fracture_width_[boundary_cells[i]], E_, nu_) - dist[i]); @@ -719,10 +770,10 @@ Fracture::solve() for (size_t i = 0; i != N; ++i) total_bnode_disp[i] += bnode_disp[i]; - + for (size_t i = 0; i != N; ++i) displacements[i] = bnode_normals_orig[i] * bnode_disp[i]; - + for (int i = 0; i != K1.size(); ++i){ //std::cout << "K1 size: " << K1.size() << std::endl; k1log << K1[i] << " "; @@ -736,18 +787,18 @@ Fracture::solve() // // @@ to facilitate debugging: identiy largest (absolute) displacement // double largest_disp = bnode_disp[0]; // //double maxK = 0; - // for (int i = 0; i != bnode_disp.size(); ++i) + // for (int i = 0; i != bnode_disp.size(); ++i) // largest_disp = abs(largest_disp) < abs(bnode_disp[i]) ? bnode_disp[i] : largest_disp; // std::cout << "max change: " << largest_disp << std::endl; // std::cout << "Max K: " << *max_element(K1.begin(), K1.end()) << std::endl; // std::cout << "Min d: " << *min_element(bnode_disp.begin(), bnode_disp.end()) << std::endl; - + bool finished = (*max_element(K1.begin(), K1.end()) <= K1max) && (*min_element(bnode_disp.begin(), bnode_disp.end()) >= -tol); - + if (finished) break; @@ -756,18 +807,18 @@ Fracture::solve() grid_stretcher_->applyBoundaryNodeDisplacements(displacements); // grid has changed its geometry, so we have to recompute discretizations - updateGridDiscretizations(); + updateGridDiscretizations(); } std::cout << "Finished! " << count << std::endl; - + // report on result }else{ OPM_THROW(std::runtime_error,"Unknowns solution method"); } } -void -Fracture::setSource() + +void Fracture::setSource() { if (rhs_pressure_.size() == 0) { size_t nc = grid_->leafGridView().size(0); @@ -804,11 +855,12 @@ Fracture::setSource() } } -double Fracture::injectionPressure() const{ - std::string control_type = prm_.get("control.type"); +double Fracture::injectionPressure() const +{ + const auto control_type = prm_.get("control.type"); if (control_type == "rate") { double bhp = 0.0; - double scale = well_source_.size(); + const double scale = well_source_.size(); // could have corrected for WI for (auto cell : well_source_) { bhp += fracture_pressure_[cell] / scale; @@ -817,97 +869,114 @@ double Fracture::injectionPressure() const{ // } else if (control_type == "bhp") { // double bhp = prm_.get("control.bhp"); } else if (control_type == "pressure") { - double bhp = prm_.get("control.pressure"); - return bhp; + return prm_.get("control.pressure"); } else if (control_type == "perf_pressure") { - double bhp = perf_pressure_; - return bhp; + return perf_pressure_; } else { OPM_THROW(std::runtime_error, "Unknowns control"); } return 0.0; } -std::vector Fracture::leakOfRate() const{ - if(leakof_.size()==0){ - std::vector rate; - return rate; +std::vector Fracture::leakOfRate() const +{ + if (leakof_.empty()) { + return {}; } + + std::vector leakofrate(grid_->leafGridView().size(0), 0.0); + ElementMapper mapper(grid_->leafGridView(), Dune::mcmgElementLayout()); - std::vector leakofrate(grid_->leafGridView().size(0),0); - for (auto& element : Dune::elements(grid_->leafGridView())) { - int eIdx = mapper.index(element); - auto geom = element.geometry(); - double dp = (fracture_pressure_[eIdx]-reservoir_pressure_[eIdx]); - double q = leakof_[eIdx]*dp; - leakofrate[eIdx] = q/geom.volume(); + for (const auto& element : Dune::elements(grid_->leafGridView())) { + const int eIdx = mapper.index(element); + const double dp = fracture_pressure_[eIdx] - reservoir_pressure_[eIdx]; + + leakofrate[eIdx] = leakof_[eIdx] * dp / element.geometry().volume(); } + return leakofrate; } -std::vector> Fracture::wellIndices() const{ - // find unique reservoir cells - if(leakof_.size() == 0){ +std::vector> Fracture::wellIndices() const +{ + if (leakof_.empty()){ // if pressure is not assembled return empty - std::vector> wellind; - return wellind; + return {}; } - std::vector res_cells = reservoir_cells_; - std::sort(res_cells.begin(),res_cells.end()); - auto last = std::unique(res_cells.begin(),res_cells.end()); - res_cells.erase(last, res_cells.end()); - std::vector q_cells(res_cells.size(),0.0); - std::vector p_cells(res_cells.size(),0.0); - double q_prev = 0; + + // find unique reservoir cells + const auto res_cells = [this]() { + auto rcells = this->reservoir_cells_; + + std::sort(rcells.begin(), rcells.end()); + rcells.erase(std::unique(rcells.begin(), rcells.end()), + rcells.end()); + + return rcells; + }(); + + std::vector q_cells(res_cells.size(), 0.0); + std::vector p_cells(res_cells.size(), 0.0); + + double q_prev = 0.0; ElementMapper mapper(grid_->leafGridView(), Dune::mcmgElementLayout()); for (auto& element : Dune::elements(grid_->leafGridView())) { - int eIdx = mapper.index(element); - double dp = (fracture_pressure_[eIdx]-reservoir_pressure_[eIdx]); - double q = leakof_[eIdx]*dp; - int res_cell = reservoir_cells_[eIdx]; - // just to search - auto it = std::find(res_cells.begin(), res_cells.end(), res_cell); - int ind_wellIdx = it - res_cells.begin(); - assert(!(it == res_cells.end())); - if( q_prev*q < 0 ){ - OPM_THROW(std::runtime_error,"Cross flow in fracture ??"); + const int eIdx = mapper.index(element); + const double dp = fracture_pressure_[eIdx] - reservoir_pressure_[eIdx]; + const double q = leakof_[eIdx] * dp; + + if (q_prev * q < 0.0) { + OPM_THROW(std::runtime_error, "Cross flow in fracture ??"); } + + // just to search + auto it = std::find(res_cells.begin(), res_cells.end(), + this->reservoir_cells_[eIdx]); + assert(it != res_cells.end()); + + const int ind_wellIdx = std::distance(res_cells.begin(), it); q_cells[ind_wellIdx] += q; p_cells[ind_wellIdx] = reservoir_pressure_[eIdx];// is set multiple times + + q_prev = q; } - std::vector> wellIndices(res_cells.size()); - double inj_press = injectionPressure(); - for(size_t i=0; i < res_cells.size(); ++i){ - double dp = inj_press - p_cells[i]; - // simplest possible approach - // assumes leakof is assumed to be calculated with reservoir cell as reference - double well_index = q_cells[i]/dp; - //assert(well_index>0); + + auto wellIndices = std::vector>{}; + wellIndices.reserve(res_cells.size()); + + const double inj_press = this->injectionPressure(); + for (auto i = 0*res_cells.size(); i < res_cells.size(); ++i) { + const double dp = inj_press - p_cells[i]; + + // Simplest possible approach. Assumes leakof is calculated with + // reservoir cell as reference. + const double well_index = q_cells[i] / dp; assert(std::isfinite(well_index)); - wellIndices[i] = {res_cells[i],well_index, origo_[2]}; - //wellIndices[i] = std::tuple({res_cells[i],well_index}); + + wellIndices.emplace_back(res_cells[i], well_index, origo_[2]); } + return wellIndices; } -void -Fracture::writePressureSystem() const{ +void Fracture::writePressureSystem() const +{ if(prm_.get("write_pressure_system")){ Dune::storeMatrixMarket(*pressure_matrix_, "pressure_matrix"); Dune::storeMatrixMarket(rhs_pressure_, "pressure_rhs"); } } -void -Fracture::writeFractureSystem() const{ +void Fracture::writeFractureSystem() const +{ if(prm_.get("write_fracture_system")){ //Dune::storeMatrixMarket(*fracture_matrix_, "fracture_matrix"); Dune::storeMatrixMarket(rhs_width_, "rhs_width"); } } -void -Fracture::solvePressure() { +void Fracture::solvePressure() +{ size_t nc = grid_->leafGridView().size(0); fracture_pressure_.resize(nc); fracture_pressure_ = 1e5; @@ -936,20 +1005,15 @@ Fracture::solvePressure() { void Fracture::solveFractureWidth() { //this->assembleFracture(); @@ should already have been taken care of - fracture_matrix_->solve(fracture_width_,rhs_width_); - double max_width = prm_.get("solver.max_width"); - double min_width = prm_.get("solver.min_width"); - for(int i=0; i < fracture_width_.size(); ++i){ - assert(std::isfinite(fracture_width_[i])); - if(fracture_width_[i]> max_width){ - std::cout << "Limit Fracture width" << std::endl; - fracture_width_[i] = max_width; - } - if(fracture_width_[i] < min_width){ - std::cout << "Remove small Fracture width" << std::endl; - fracture_width_[i] = min_width; - } - assert(std::isfinite(fracture_width_[i])); + fracture_matrix_->solve(fracture_width_, rhs_width_); + + const double max_width = prm_.get("solver.max_width"); + const double min_width = prm_.get("solver.min_width"); + + for (auto& width : this->fracture_width_) { + width = std::clamp(width[0], min_width, max_width); + + assert (std::isfinite(width)); } } @@ -959,12 +1023,11 @@ void Fracture::initFractureStates() this->initFracturePressureFromReservoir(); } -void -Fracture::initFractureWidth() +void Fracture::initFractureWidth() { - size_t nc = grid_->leafGridView().size(0); - fracture_width_.resize(nc); + fracture_width_.resize(grid_->leafGridView().size(0)); fracture_width_ = prm_.get("config.initial_fracture_width"); + // ElementMapper elemMapper(grid_->leafGridView(), Dune::mcmgElementLayout()); // for (auto& element : elements(grid_->leafGridView())) { // const auto elemIdx = elemMapper.index(element); @@ -974,67 +1037,65 @@ Fracture::initFractureWidth() // fracture_width_[elemIdx] *= dist_origo; // } } -void -Fracture::initFracturePressureFromReservoir(){ - size_t nc = reservoir_cells_.size(); - fracture_pressure_.resize(nc); - for(size_t i=0; i < nc; ++i){ - fracture_pressure_[i] = reservoir_pressure_[i]; - } + +void Fracture::initFracturePressureFromReservoir() +{ + const auto nc = reservoir_cells_.size(); + + this->fracture_pressure_.resize(nc); + std::copy_n(this->reservoir_pressure_.begin(), nc, + this->fracture_pressure_.begin()); } -void -Fracture::initPressureMatrix() + +void Fracture::initPressureMatrix() { - // size_t num_columns = 0; - // index of the neighbour - // - double fWI = prm_.get("fractureWI"); - for(int cell : well_source_){ - perfinj_.push_back({cell,fWI}); + { + const double fWI = prm_.get("fractureWI"); + + for (const auto& cell : well_source_) { + perfinj_.emplace_back(cell, fWI); + } } - size_t nc = grid_->leafGridView().size(0); - leakof_.resize(nc,0.0); + const auto nc = grid_->leafGridView().size(0); + + leakof_.resize(nc, 0.0); + ElementMapper mapper(grid_->leafGridView(), Dune::mcmgElementLayout()); - for (auto& element : Dune::elements(grid_->leafGridView())) { - int eIdx = mapper.index(element); - auto geom = element.geometry(); + for (const auto& element : Dune::elements(grid_->leafGridView())) { + const int eIdx = mapper.index(element); + const auto eCenter = element.geometry().center(); + // iterator over all intersections - for (auto& is : Dune::intersections(grid_->leafGridView(),element)) { - - if (!is.boundary()) { - auto eCenter = geom.center(); - int nIdx = mapper.index(is.outside()); - if (eIdx < nIdx) { - // calculate distance between the midpoints - auto nCenter = is.outside().geometry().center(); - auto isCenter = is.geometry().center(); - auto d_inside = eCenter - isCenter; - auto d_outside = nCenter - isCenter; - // probably should use projected distenace - auto igeom = is.geometry(); - double area = igeom.volume(); - double h1 = area/d_inside.two_norm(); - double h2 = area/d_outside.two_norm(); - - { - Htrans matel(nIdx, eIdx, h1, h2); - htrans_.push_back(matel); - } - } + for (const auto& is : Dune::intersections(grid_->leafGridView(), element)) { + if (is.boundary()) { + continue; } + + const int nIdx = mapper.index(is.outside()); + if (eIdx >= nIdx) { + continue; + } + + // calculate distance between the midpoints + const auto nCenter = is.outside().geometry().center(); + const auto isCenter = is.geometry().center(); + const auto d_inside = eCenter - isCenter; + const auto d_outside = nCenter - isCenter; + + // probably should use projected distenace + const double area = is.geometry().volume(); + const double h1 = area / d_inside.two_norm(); + const double h2 = area / d_outside.two_norm(); + + this->htrans_.emplace_back(nIdx, eIdx, h1, h2); } - { - //auto normal = ddm::normalOfElement(element); - //auto permmatrix = reservoir_perm_[eIdx]; - //auto pn = permmatrix.mv(normal); - //double permval = pn.dot(normal); - - // keap reservoir perm as n'K'n - double value = reservoir_mobility_[eIdx]*(reservoir_perm_[eIdx] * geom.volume()) ; - value /= reservoir_dist_[eIdx]; - leakof_[eIdx] = value; - } + + // Keep reservoir perm as n'K'n. + leakof_[eIdx] = reservoir_mobility_[eIdx] + * reservoir_perm_[eIdx] + * element.geometry().volume() + / reservoir_dist_[eIdx]; } // not need if build mode is implicit @@ -1049,10 +1110,10 @@ Fracture::initPressureMatrix() //matrix.setImplicitBuildModeParameters(3 * 6 - 3 - 2, 0.4); //size_t nc = grid_->leafGridView().size(0); //matrix.setSize(nc, nc); - for (auto matel : htrans_) { - size_t i = std::get<0>(matel); - size_t j = std::get<1>(matel); - double zero_entry = 0.0; //1e-11; + for (const auto& matel : htrans_) { + const auto i = std::get<0>(matel); + const auto j = std::get<1>(matel); + const double zero_entry = 0.0; //1e-11; matrix.entry(i, j) = zero_entry; //0; matrix.entry(j, i) = zero_entry; //0; matrix.entry(j, j) = zero_entry; //0; @@ -1061,25 +1122,31 @@ Fracture::initPressureMatrix() matrix.compress(); //} } -void -Fracture::assemblePressure() + +void Fracture::assemblePressure() { - auto& matrix = *pressure_matrix_; + const auto control_type = prm_.get_child("control").get("type"); + const auto isPressCtrl = (control_type == "pressure") || + (control_type == "perf_pressure"); + + if ((control_type != "rate") && !isPressCtrl) { + OPM_THROW(std::runtime_error, "Unknown control of injection into Fracture"); + } + + auto& matrix = *this->pressure_matrix_; matrix = 0.0; - double mobility=1e4; //1e4; // @@ 1.0 - - for (auto matel : htrans_) { - size_t i = std::get<0>(matel); - size_t j = std::get<1>(matel); - double t1 = std::get<2>(matel); - double t2 = std::get<3>(matel); - double h1 = fracture_width_[i]; - double h2 = fracture_width_[j]; + + const double mobility = 1.0e4; //1e4; // @@ 1.0 + + for (const auto& [i, j, t1, t2] : this->htrans_) { + const double h1 = fracture_width_[i]; + const double h2 = fracture_width_[j]; + // harmonic mean of surface flow double value = 12. / (h1 * h1 * h1 * t1) + 12. / (h2 * h2 * h2 * t2); - value = 1 / value; value *= mobility; + // matrix.entry(i, j) -= value; // matrix.entry(j, i) -= value; // // @@ -1092,27 +1159,23 @@ Fracture::assemblePressure() matrix[i][i] += value; matrix[j][j] += value; } - auto control = prm_.get_child("control"); - std::string control_type = control.get("type"); - for (size_t i = 0; i < leakof_.size(); ++i) { - // matrix.entry(i, i) += leakof_[i]; - matrix[i][i] += leakof_[i]; + + for (auto i = 0*leakof_.size(); i < leakof_.size(); ++i) { + // matrix.entry(i, i) += leakof_[i]; + matrix[i][i] += leakof_[i]; } - if (control_type == "rate") { - // no extra tings in matrix - } else if (control_type == "pressure" || "perf_pressure") { - for (const auto& perfinj : perfinj_) { - int cell = std::get<0>(perfinj); - double value = std::get<1>(perfinj); + + if (isPressCtrl) { + for (const auto& [cell, value] : perfinj_) { matrix[cell][cell] += value; } - }else{ - OPM_THROW(std::runtime_error,"Unknown control of injection into Fracture"); } } -void Fracture::assembleFracture(){ - size_t nc = grid_->leafGridView().size(0); +void Fracture::assembleFracture() +{ + const size_t nc = grid_->leafGridView().size(0); + //rhs_width_.resize(nc); rhs_width_ = fracture_pressure_; ElementMapper mapper(grid_->leafGridView(), Dune::mcmgElementLayout()); @@ -1126,21 +1189,22 @@ void Fracture::assembleFracture(){ // using compressible stress ?? rhs_width_[idx] -= ddm::tractionSymTensor(reservoir_stress_[idx],cell_normals_[idx]); - + // @@ not entirely accurate, but will avoid unphysical negative normal displacements if (rhs_width_[idx] < 0) rhs_width_[idx] = 0; } - - } + // Do we need to ad thermal "forces" if(!fracture_matrix_){ fracture_matrix_ = std::make_unique(); } + fracture_matrix_->resize(nc,nc); *fracture_matrix_ = 0.0; + ddm::assembleMatrix(*fracture_matrix_,E_, nu_,*grid_); } @@ -1184,11 +1248,11 @@ void Fracture::printMechMatrix() const // debug purposes // } - template void -Fracture::updateReservoirCells(const external::cvf::ref& cellSearchTree, - const Dune::CpGrid& grid3D); +Fracture::updateReservoirCells(const external::cvf::ref& cellSearchTree, + const Dune::CpGrid& grid3D); template void -Fracture::updateReservoirCells>(const external::cvf::ref& cellSearchTree, - const Dune::PolyhedralGrid<3,3,double>& grid3D); +Fracture::updateReservoirCells(const external::cvf::ref& cellSearchTree, + const Dune::PolyhedralGrid<3,3,double>& grid3D); + } // namespace Opm diff --git a/opm/geomech/Fracture.hpp b/opm/geomech/Fracture.hpp index 92d2a39..afeec53 100644 --- a/opm/geomech/Fracture.hpp +++ b/opm/geomech/Fracture.hpp @@ -176,12 +176,14 @@ class Fracture void setFractureGrid(std::unique_ptr& gptr); // a hack to allow use of another grid std::vector> wellIndices() const; WellInfo& wellInfo(){return wellinfo_;} + const WellInfo& wellInfo() const {return wellinfo_;} std::vector leakOfRate() const; double injectionPressure() const; void setPerfPressure(double perfpressure){perf_pressure_ = perfpressure;} - Dune::FieldVector stress(Dune::FieldVector obs) const; - Dune::FieldVector strain(Dune::FieldVector obs) const; - Dune::FieldVector disp(Dune::FieldVector obs) const; + Dune::FieldVector stress(const Dune::FieldVector& obs) const; + Dune::FieldVector strain(const Dune::FieldVector& obs) const; + Dune::FieldVector disp(const Dune::FieldVector& obs) const; + private: Dune::BlockVector> all_slips() const; void resetWriters(); diff --git a/opm/geomech/FractureModel.cpp b/opm/geomech/FractureModel.cpp index 2fb14cb..3481164 100644 --- a/opm/geomech/FractureModel.cpp +++ b/opm/geomech/FractureModel.cpp @@ -1,206 +1,257 @@ #include "config.h" + #include -#include #include + +#include + #include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + namespace Opm{ - void FractureModel::addWell(std::string name, - const std::vector& points, - const std::vector& segments, - const std::vector& res_cells){ - std::string outputdir = prm_.get("outputdir"); - std::string casename = prm_.get("casename"); - wells_.push_back(FractureWell(outputdir,casename, name,points,segments, res_cells)); - // add witth no fractures - well_fractures_.push_back(std::vector()); + void FractureModel::addWell(const std::string& name, + const std::vector& points, + const std::vector& segments, + const std::vector& res_cells) + { + const std::string outputdir = prm_.get("outputdir"); + const std::string casename = prm_.get("casename"); + + wells_.emplace_back(outputdir, casename, name, points, segments, res_cells); + + // add with no fractures + well_fractures_.emplace_back(); } - void - FractureModel::addFractures() + void FractureModel::addFractures() { - auto config = prm_.get_child("config"); - std::string type = config.get("type"); - for (size_t i = 0; i < wells_.size(); ++i) { - const FractureWell& well = wells_[i]; - auto& grid = well.grid(); - using ElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper; + const auto config = prm_.get_child("config"); + const std::string type = config.get("type"); + + using ElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper + ; + + auto wellFrac = this->well_fractures_.begin(); + for (const auto& well : this->wells_) { + const auto& grid = well.grid(); + ElementMapper mapper(grid.leafGridView(), Dune::mcmgElementLayout()); for (const auto& elem : elements(grid.leafGridView())) { - int eIdx = mapper.index(elem); - auto geo = elem.geometry(); + const int eIdx = mapper.index(elem); + const auto geo = elem.geometry(); assert(geo.corners() == 2); + // auto origo = geo.center(); Dune::FieldVector origo, normal; int perf = eIdx; - int well_cell = well.reservoirCell(eIdx); if (type == "perp_well") { - origo = geo.corner(1); // assume this is cell center - normal = geo.corner(1) - geo.corner(0); - // perf = well_fractures_[i].size(); - } else if (type == "well_seed") { - if( config.get("well") == well.name()){ - std::cout << "Fracure added for well " << well.name() << std::endl; - //std::vector cell_ijk = config.get< std::vector > ("cell_ijk"); - int cell = config.get< int> ("cell"); - if(well.reservoirCell(eIdx) == cell){ - origo = geo.corner(1); - auto config_bst = config.getBoostParamPtr(); - //double tmp = config_bst->get ("normal",1); - //for (auto i : as_vector(pt, "a")) std::cout << i << ' '; - std::vector tmp_normal = as_vector(*config_bst,"normal"); - assert(tmp_normal.size() == 3); // wrong use of tmpmal. - for(int i=0; i < 3; ++i){ - normal[i] = tmp_normal[i]; - } - }else{ - continue; - } - }else{ - continue; - } - } else if (type == "tensile_fracture") { - // https://link.springer.com/article/10.1007/s40948-023-00694-1 - double fractureThoughness = 1.0e6; // reservoir_fractureThoughness_[eIdx]; // ca 1.0 MPa m^1/2 - double tensilestrength = 5e6; // reservoir_tensilestrength_[eIdx]; // 5 MPa - double criticallength = (fractureThoughness / tensilestrength); // ca (1/5)^2 5 mm. - criticallength *= criticallength; - auto stressmat = ddm::symTensor2Matrix(well.reservoirStress(eIdx)); - Dune::FieldMatrix eigenVectors; - Dune::FieldVector eigenValues; - Dune::FMatrixHelp::eigenValuesVectors(stressmat, eigenValues, eigenVectors); - int min_dir = -1; - int max_dir = -1; - double min_eig = 1e99; - double max_eig = -1e99; - for (int i = 0; i < 3; ++i) { - if (eigenValues[i] < min_eig) { - min_dir = i; - min_eig = eigenValues[i]; - } - if (eigenValues[i] > max_eig) { - max_dir = i; - max_eig = eigenValues[i]; - } - } - normal = eigenVectors[min_dir]; - // take midpoint - origo = geo.corner(1); //-geo.corner(0); - // origo /= 2.0; - // expression for size; - } else { - OPM_THROW(std::runtime_error, "Invalid fracture type"); + origo = geo.corner(1); // assume this is cell center + normal = geo.corner(1) - geo.corner(0); + // perf = well_fractures_[i].size(); + } + else if (type == "well_seed") { + if (config.get("well") == well.name()) { + std::cout << "Fracure added for well " << well.name() << std::endl; + //std::vector cell_ijk = config.get< std::vector > ("cell_ijk"); + int cell = config.get< int> ("cell"); + if (well.reservoirCell(eIdx) == cell) { + origo = geo.corner(1); + auto config_bst = config.getBoostParamPtr(); + //double tmp = config_bst->get ("normal",1); + //for (auto i : as_vector(pt, "a")) std::cout << i << ' '; + std::vector tmp_normal = as_vector(*config_bst,"normal"); + assert(tmp_normal.size() == 3); // wrong use of tmpmal. + for (int i=0; i < 3; ++i) { + normal[i] = tmp_normal[i]; + } + } + else { + continue; + } + } + else { + continue; + } } + else if (type == "tensile_fracture") { + // https://link.springer.com/article/10.1007/s40948-023-00694-1 + const double fractureThoughness = 1.0e6; // reservoir_fractureThoughness_[eIdx]; // ca 1.0 MPa m^1/2 + const double tensilestrength = 5e6; // reservoir_tensilestrength_[eIdx]; // 5 MPa + + double criticallength = (fractureThoughness / tensilestrength); // ca (1/5)^2 5 mm. + criticallength *= criticallength; + + Dune::FieldMatrix eigenVectors; + Dune::FieldVector eigenValues; + auto stressmat = ddm::symTensor2Matrix(well.reservoirStress(eIdx)); + Dune::FMatrixHelp::eigenValuesVectors(stressmat, eigenValues, eigenVectors); + + int min_dir = -1; + int max_dir = -1; + double min_eig = 1e99; + double max_eig = -1e99; + for (int i = 0; i < 3; ++i) { + if (eigenValues[i] < min_eig) { + min_dir = i; + min_eig = eigenValues[i]; + } + if (eigenValues[i] > max_eig) { + max_dir = i; + max_eig = eigenValues[i]; + } + } + + normal = eigenVectors[min_dir]; + // take midpoint + origo = geo.corner(1); //-geo.corner(0); + // origo /= 2.0; + // expression for size; + } + else { + OPM_THROW(std::runtime_error, "Invalid fracture type"); + } - Fracture fracture; - fracture.init(well.name(), perf, well_cell, origo, normal, prm_); - well_fractures_[i].push_back(std::move(fracture)); + wellFrac->emplace_back() + .init(well.name(), perf, well_cell, origo, normal, prm_); } + + ++wellFrac; } } - void FractureModel::initFractureStates(){ - for(size_t i=0; i < wells_.size(); ++i){ - std::vector& fractures = well_fractures_[i]; - for(size_t j=0; j < fractures.size(); ++j){ - fractures[j].initFractureStates(); + void FractureModel::initFractureStates() + { + for (auto& fractures : this->well_fractures_) { + for (auto& fracture : fractures) { + fracture.initFractureStates(); } } } + // probably this should be collected in one loop - Dune::FieldVector FractureModel::stress(Dune::FieldVector obs) const{ - Dune::FieldVector stress; - for(size_t i=0; i < wells_.size(); ++i){ - const std::vector& fractures = well_fractures_[i]; - for(size_t j=0; j < fractures.size(); ++j){ - Dune::FieldVector tmp = fractures[j].stress(obs); - stress += tmp; + Dune::FieldVector + FractureModel::stress(const Dune::FieldVector& obs) const + { + Dune::FieldVector stress{}; + + for (const auto& fractures : this->well_fractures_) { + for (const auto& fracture : fractures) { + stress += fracture.stress(obs); } } + return stress; } - Dune::FieldVector FractureModel::strain(Dune::FieldVector obs) const{ - Dune::FieldVector strain; - for(size_t i=0; i < wells_.size(); ++i){ - const std::vector& fractures = well_fractures_[i]; - for(size_t j=0; j < fractures.size(); ++j){ - Dune::FieldVector tmp = fractures[j].strain(obs); - strain += tmp; + + Dune::FieldVector + FractureModel::strain(const Dune::FieldVector& obs) const + { + Dune::FieldVector strain{}; + + for (const auto& fractures : this->well_fractures_) { + for (const auto& fracture : fractures) { + strain += fracture.strain(obs); } } + return strain; } - Dune::FieldVector FractureModel::disp(Dune::FieldVector obs) const{ - Dune::FieldVector disp; - for(size_t i=0; i < wells_.size(); ++i){ - const std::vector& fractures = well_fractures_[i]; - for(size_t j=0; j < fractures.size(); ++j){ - Dune::FieldVector tmp = fractures[j].disp(obs); - disp += tmp; + + Dune::FieldVector + FractureModel::disp(const Dune::FieldVector& obs) const + { + Dune::FieldVector disp{}; + + for (const auto& fractures : this->well_fractures_) { + for (const auto& fracture : fractures) { + disp += fracture.disp(obs); } } - return disp; + return disp; } - void FractureModel::write(int reportStep) const{ - for(size_t i=0; i < wells_.size(); ++i){ - wells_[i].write(); - const std::vector& fractures = well_fractures_[i]; - for(size_t j=0; j < fractures.size(); ++j){ - fractures[j].write(reportStep); + void FractureModel::write(int reportStep) const + { + for (auto i = 0*this->wells_.size(); i < this->wells_.size(); ++i) { + this->wells_[i].write(); + + for (const auto& fracture : this->well_fractures_[i]) { + fracture.write(reportStep); } } } - void FractureModel::writemulti(double time) const{ - for(size_t i=0; i < wells_.size(); ++i){ - wells_[i].writemulti(time); - const std::vector& fractures = well_fractures_[i]; - for(size_t j=0; j < fractures.size(); ++j){ - fractures[j].writemulti(time); + + void FractureModel::writemulti(double time) const + { + for (auto i = 0*this->wells_.size(); i < this->wells_.size(); ++i) { + this->wells_[i].writemulti(time); + + for (const auto& fracture : this->well_fractures_[i]) { + fracture.writemulti(time); } } } - void FractureModel::solve() { - for(size_t i=0; i < wells_.size(); ++i){ - std::vector& fractures = well_fractures_[i]; - for(size_t j=0; j < fractures.size(); ++j){ - fractures[j].solve(); + + void FractureModel::solve() + { + for (auto& fractures : this->well_fractures_) { + for (auto& fracture : fractures) { + fracture.solve(); } } } - void FractureModel::updateReservoirProperties() { - for(size_t i=0; i < wells_.size(); ++i){ - std::vector& fractures = well_fractures_[i]; - for(size_t j=0; j < fractures.size(); ++j){ - fractures[j].updateReservoirProperties(); + + void FractureModel::updateReservoirProperties() + { + for (auto& fractures : this->well_fractures_) { + for (auto& fracture : fractures) { + fracture.updateReservoirProperties(); } } } + std::vector> - FractureModel::getExtraWellIndices(std::string wellname) const{ - // for now just do a search - bool addconnections = prm_.get("addconnections"); - if (addconnections) { - for (size_t i = 0; i < wells_.size(); ++i) { - if (wells_[i].name() == wellname) { - // collect all from a well - std::vector> wellindices; - for (const auto& frac : well_fractures_[i]) { - auto perfs = frac.wellIndices(); - wellindices.insert(wellindices.begin(),perfs.begin(), perfs.end()); - } - return wellindices; - } - } - std::string message = "Now fractures on this well found"; - message += wellname; - OPM_THROW(std::runtime_error, message.c_str()); + FractureModel::getExtraWellIndices(const std::string& wellname) const + { + auto wellindices = std::vector>{}; + + if (! prm_.get("addconnections")) { + return wellindices; } - std::vector> empty; - return empty; + + // For now just do a search + auto pos = std::find_if(this->wells_.begin(), this->wells_.end(), + [&wellname](const auto& well) + { return well.name() == wellname; }); + + if (pos == this->wells_.end()) { + OPM_THROW(std::runtime_error, "No fractures found for well " + wellname); + } + + const auto i = std::distance(this->wells_.begin(), pos); + for (const auto& frac : well_fractures_[i]) { + auto perfs = frac.wellIndices(); + + wellindices.insert(wellindices.begin(), + std::make_move_iterator(perfs.begin()), + std::make_move_iterator(perfs.end())); + } + + return wellindices; } } diff --git a/opm/geomech/FractureModel.hpp b/opm/geomech/FractureModel.hpp index 9544525..35946e2 100644 --- a/opm/geomech/FractureModel.hpp +++ b/opm/geomech/FractureModel.hpp @@ -2,24 +2,32 @@ #include "Fracture.hpp" #include "FractureWell.hpp" #include "GeometryHelpers.hpp" + #include #include #include + #include #include #include + #include + #include #include #include + #include + using ptree = boost::property_tree::ptree; template std::vector as_vector(ptree const& pt, ptree::key_type const& key) { std::vector r; - for (auto& item : pt.get_child(key)) + for (const auto& item : pt.get_child(key)) { r.push_back(item.second.get_value()); + } + return r; } @@ -31,8 +39,9 @@ std::vector as_vector(ptree const& pt, ptree::key_type const& key) // r.push_back(item.second.get_value()); // return r; // } -namespace Opm{ -class FractureModel{ +namespace Opm { +class FractureModel +{ //using CartesianIndexMapper = Dune::CartesianIndexMapper; public: using Point3D = Dune::FieldVector; @@ -40,8 +49,8 @@ class FractureModel{ template FractureModel(const Grid& grid, const std::vector& wells, - const Opm::PropertyTree& - ); + const Opm::PropertyTree& param); + void addFractures(); template @@ -75,18 +84,24 @@ class FractureModel{ } } } - void addWell(std::string name, const std::vector& points,const std::vector>& segments,const std::vector& reservoir_cells ); + + void addWell(const std::string& name, + const std::vector& points, + const std::vector>& segments, + const std::vector& reservoir_cells); + void write(int ReportStep = -1) const; void writemulti(double time) const; void solve(); void updateReservoirProperties(); void initFractureStates(); + template void initReservoirProperties(const Simulator& simulator) { - for (size_t i=0; i < wells_.size(); ++i) { - for (auto& fracture : well_fractures_[i]){ - fracture.initReservoirProperties(simulator); + for (auto& fractures : this->well_fractures_) { + for (auto& fracture : fractures) { + fracture.initReservoirProperties(simulator); } } } @@ -96,21 +111,20 @@ class FractureModel{ { for (size_t i=0; i < wells_.size(); ++i) { for (auto& fracture : well_fractures_[i]){ - fracture.updateReservoirProperties(simulator); + fracture.updateReservoirProperties(simulator); // set well properties - WellInfo wellinfo = fracture.wellInfo(); + const auto& wellinfo = fracture.wellInfo(); // need to be double checked how to assosiate correct perforation/segment - int perf_index_frac = wellinfo.perf; - int cell_index_frac = wellinfo.well_cell; - std::size_t well_index = simulator.problem().wellModel().wellState().index(wellinfo.name).value(); - const auto& wellstate = simulator.problem().wellModel().wellState().well(well_index); + const int perf_index_frac = wellinfo.perf; + const int cell_index_frac = wellinfo.well_cell; + const auto& wellstate = simulator.problem().wellModel().wellState().well(wellinfo.name); const auto& perf_data = wellstate.perf_data; auto it = std::find(perf_data.cell_index.begin(), perf_data.cell_index.end(), cell_index_frac); - int perf_index = it- perf_data.cell_index.begin(); + const int perf_index = it- perf_data.cell_index.begin(); std::cout << "Perf index flow " << perf_index << " fracture " << perf_index_frac << std::endl; - double perf_pressure = perf_data.pressure[perf_index]; + const double perf_pressure = perf_data.pressure[perf_index]; fracture.setPerfPressure(perf_pressure); wells_[i].setPerfPressure(perf_index_frac, perf_pressure); //NB do we need some rates? need to be summed over "peforations of the fractures" @@ -118,22 +132,29 @@ class FractureModel{ } } - template - void updateReservoirWellProperties(const Simulator& simulator) { - for(auto& well : wells_){ - well.updateReservoirProperties(simulator); + template + void updateReservoirWellProperties(const Simulator& simulator) + { + for (auto& well : wells_) { + well.updateReservoirProperties(simulator); } } - std::vector> getExtraWellIndices(std::string wellname) const; + + std::vector> + getExtraWellIndices(const std::string& wellname) const; + bool addPertsToSchedule(){return prm_.get("addperfs_to_schedule");} + // probably this should be collected in one loop since all do full loop over fracture ++ well - Dune::FieldVector stress(Dune::FieldVector obs) const; - Dune::FieldVector strain(Dune::FieldVector obs) const; - Dune::FieldVector disp(Dune::FieldVector obs) const; + Dune::FieldVector stress(const Dune::FieldVector& obs) const; + Dune::FieldVector strain(const Dune::FieldVector& obs) const; + Dune::FieldVector disp(const Dune::FieldVector& obs) const; + private: std::vector wells_; std::vector> well_fractures_; Opm::PropertyTree prm_; }; } + #include "FractureModel_impl.hpp" diff --git a/opm/geomech/FractureModel_impl.hpp b/opm/geomech/FractureModel_impl.hpp index c7da22d..7ddbc98 100644 --- a/opm/geomech/FractureModel_impl.hpp +++ b/opm/geomech/FractureModel_impl.hpp @@ -3,83 +3,81 @@ namespace Opm{ template FractureModel::FractureModel(const Grid& grid, - const std::vector& wells, - const Opm::PropertyTree& param) - : - prm_(param) + const std::vector& wells, + const Opm::PropertyTree& param) + : prm_(param) { - using CartesianIndexMapper = Dune::CartesianIndexMapper; - CartesianIndexMapper cartmapper(grid); - const auto& cartdim = cartmapper.cartesianDimensions(); - std::vector cart(cartdim[0]*cartdim[1]*cartdim[2], -1); - using ElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper; - ElementMapper mapper(grid.leafGridView(), Dune::mcmgElementLayout()); - for (const auto& elem : elements(grid.leafGridView())) { - int eIdx = mapper.index(elem); - cart[ cartmapper.cartesianIndex( eIdx ) ] = eIdx; - } - const auto& config = prm_.get_child("config"); - std::string type = config.get("type"); - if( type == "well_seed"){ - auto config_bst = prm_.getBoostParamPtr(); - std::vector cell_ijk = as_vector(*config_bst, "config.cell_ijk"); - std::array ijk; - assert(cell_ijk.size() == 3); - for(int i=0; i < 3; ++i){ - // map to 0 based - ijk[i] = cell_ijk[i]-1; - } - int cartIdx = ijk[0]+(cartdim[0])*ijk[1]+(cartdim[0]*cartdim[1])*ijk[2];//NB assumes ijk ordering - int reservoir_cell = cart[cartIdx]; - assert(cartIdx>=0); - prm_.put("config.cell",reservoir_cell); - } - - + using CartesianIndexMapper = Dune::CartesianIndexMapper; + CartesianIndexMapper cartmapper(grid); + const auto& cartdim = cartmapper.cartesianDimensions(); + std::vector cart(cartdim[0]*cartdim[1]*cartdim[2], -1); + using ElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper; + ElementMapper mapper(grid.leafGridView(), Dune::mcmgElementLayout()); + for (const auto& elem : elements(grid.leafGridView())) { + const int eIdx = mapper.index(elem); + cart[ cartmapper.cartesianIndex(eIdx) ] = eIdx; + } + const auto& config = prm_.get_child("config"); + std::string type = config.get("type"); + if (type == "well_seed"){ + auto config_bst = prm_.getBoostParamPtr(); + std::vector cell_ijk = as_vector(*config_bst, "config.cell_ijk"); + std::array ijk; + assert(cell_ijk.size() == 3); + for(int i=0; i < 3; ++i){ + // map to 0 based + ijk[i] = cell_ijk[i]-1; + } + int cartIdx = ijk[0] + cartdim[0]*(ijk[1] + cartdim[1]*ijk[2]);//NB assumes ijk ordering + int reservoir_cell = cart[cartIdx]; + assert(cartIdx>=0); + prm_.put("config.cell",reservoir_cell); + } + //using CartesianIndexMapper = Dune::CartesianIndexMapper; //CartesianIndexMapper cartmapper(grid); //const std::array //auto cartdims = cartmapper.cartesianDimensions(); GeometryHelper geomhelp(grid); - // NB: need to be carefull in parallel + // NB: need to be carefull in parallel for(int wellIdx=0; wellIdx < wells.size(); ++wellIdx){ std::vector vertices; std::vector wellsegment; auto well = wells[wellIdx]; if(true){//well.isInjector()){// only look at injectors - std::vector cells; - // should be done property with welltraj - for(const auto& connection : well.getConnections()){ - int global_index = connection.global_index(); - int cell_index = geomhelp.compressedIndex(global_index); - Point3D vertex = geomhelp.centroid(cell_index); - vertices.push_back(vertex); - cells.push_back(cell_index); - } - Point3D refpoint(vertices[0]); - if(well.hasRefDepth()){ - if(std::abs(well.getRefDepth()-refpoint[2])< 10){ - refpoint[2] = well.getRefDepth()-10;//avoid zero + std::vector cells; + // should be done property with welltraj + for(const auto& connection : well.getConnections()){ + int global_index = connection.global_index(); + int cell_index = geomhelp.compressedIndex(global_index); + Point3D vertex = geomhelp.centroid(cell_index); + vertices.push_back(vertex); + cells.push_back(cell_index); + } + Point3D refpoint(vertices[0]); + if(well.hasRefDepth()){ + if(std::abs(well.getRefDepth()-refpoint[2])< 10){ + refpoint[2] = well.getRefDepth()-10;//avoid zero + }else{ + refpoint[2] = well.getRefDepth(); + } }else{ - refpoint[2] = well.getRefDepth(); + refpoint[2] -= 10; } - }else{ - refpoint[2] -= 10; - } - vertices.insert(vertices.begin(),refpoint); - std::vector segments; - //NB unsinged so grid can not be to big - for(unsigned i=0; i < unsigned(cells.size()); ++i){ - segments.push_back(Segment({i,i+1})); - } - // NB should add gri cells - this->addWell(well.name(),vertices, segments, cells); + vertices.insert(vertices.begin(),refpoint); + std::vector segments; + //NB unsinged so grid can not be to big + for(unsigned i=0; i < unsigned(cells.size()); ++i){ + segments.push_back(Segment({i,i+1})); + } + // NB should add gri cells + this->addWell(well.name(),vertices, segments, cells); } } } template void FractureModel::addDefaultsWells(const Grid& grid, const Opm::EclipseGrid& eclgrid){ - this->addFractures(); - this->updateFractureReservoirCells(grid,eclgrid); + this->addFractures(); + this->updateFractureReservoirCells(grid,eclgrid); } } diff --git a/opm/geomech/FractureWell.cpp b/opm/geomech/FractureWell.cpp index f0c7c6c..627eb2c 100644 --- a/opm/geomech/FractureWell.cpp +++ b/opm/geomech/FractureWell.cpp @@ -1,42 +1,58 @@ #include "config.h" + #include + #include -namespace Opm{ + +namespace Opm { void FractureWell::init(const std::vector& points, - const std::vector& segments){ - Dune::GridFactory factory; - for(size_t i=0; i < points.size(); ++i){ - factory.insertVertex(points[i]); + const std::vector& segments) + { + Dune::GridFactory factory{}; + + for (const auto& vertex : points) { + factory.insertVertex(vertex); } - for (size_t i = 0; i < segments.size(); ++i) { + + for (const auto& segment : segments) { std::vector seg(2,0); - seg[0] = segments[i][0];seg[1] = segments[i][1]; + seg[0] = segment[0]; + seg[1] = segment[1]; + factory.insertElement(Dune::GeometryTypes::line, seg); // std::move(mapping)); } + grid_ = factory.createGrid(); + this->resetWriters(); } - void FractureWell::write() const{ - std::string filename = outputprefix_ + "/" + this->name(); - vtkwriter_->write(filename.c_str()); + void FractureWell::write() const + { + const std::string filename = outputprefix_ + "/" + this->name(); + vtkwriter_->write(filename); } - void FractureWell::resetWriters(){ + void FractureWell::resetWriters() + { // nead to be reseat if grid is changed ?? - vtkwriter_ = std::make_unique>(grid_->leafGridView(), Dune::VTK::nonconforming); - std::string outputdir = outputprefix_; - std::string simName = casename_ + "_" + this->name(); - std::string multiFileName = ""; - vtkmultiwriter_ = std::make_unique< Opm::VtkMultiWriter >(/*async*/ false, - grid_->leafGridView(), - outputdir, - simName, - multiFileName - ); + vtkwriter_ = std::make_unique> + (grid_->leafGridView(), Dune::VTK::nonconforming); + + const std::string outputdir = outputprefix_; + const std::string simName = casename_ + "_" + this->name(); + const std::string multiFileName = ""; + + vtkmultiwriter_ = std::make_unique> + (/*async*/ false, + grid_->leafGridView(), + outputdir, + simName, + multiFileName); } - void FractureWell::writemulti(double time) const{ + void FractureWell::writemulti(const double time) const + { std::vector reservoir_pressure = reservoir_pressure_; std::vector reservoir_temperature = reservoir_temperature_; //std::vector reservoir_cells = reservoir_cells_; @@ -54,21 +70,25 @@ namespace Opm{ vtkmultiwriter_->endWrite(); } - FractureWell::FractureWell(std::string outputprefix, - std::string casename, - std::string name, - const std::vector& points, - const std::vector& segments, - const std::vector& res_cells){ - outputprefix_ = outputprefix; - casename_ = casename; + FractureWell::FractureWell(const std::string& outputprefix, + const std::string& casename, + const std::string& name, + const std::vector& points, + const std::vector& segments, + const std::vector& res_cells) + { + outputprefix_ = outputprefix; + casename_ = casename; name_ = name; + this->init(points, segments); + vtkwriter_ = std::make_unique> (grid_->leafGridView(), Dune::VTK::nonconforming); + reservoir_cells_ = res_cells; reservoir_stress_.resize(res_cells.size()); - std::fill(reservoir_stress_.begin(),reservoir_stress_.end(),100e5);// random initialization + std::fill(reservoir_stress_.begin(), reservoir_stress_.end(), 100.0e5); // random initialization } } diff --git a/opm/geomech/FractureWell.hpp b/opm/geomech/FractureWell.hpp index 606bb8f..aedf8e3 100644 --- a/opm/geomech/FractureWell.hpp +++ b/opm/geomech/FractureWell.hpp @@ -7,9 +7,9 @@ class FractureWell using Segment = std::array; public: using Grid = Dune::FoamGrid<1,3>; - FractureWell(std::string outputprefix, - std::string casename, - std::string name, + FractureWell(const std::string& outputprefix, + const std::string& casename, + const std::string& name, const std::vector& points, const std::vector& segments, const std::vector& res_cells); @@ -43,12 +43,18 @@ class FractureWell } void setPerfPressure(int perf_index, double pressure){perf_pressure_[perf_index] = pressure;} const Grid& grid() const{return *grid_;} - std::string name() const{return name_;}; + const std::string& name() const{return name_;}; void write() const; void writemulti(double time) const; void resetWriters(); - int reservoirCell(int wellcell) const {return reservoir_cells_[wellcell];}; - Dune::FieldVector reservoirStress(int wellcell) const{return reservoir_stress_[wellcell];}; + int reservoirCell(int wellcell) const {return reservoir_cells_[wellcell];} + + const Dune::FieldVector& + reservoirStress(int wellcell) const + { + return reservoir_stress_[wellcell]; + } + private: std::string outputprefix_; std::string casename_; diff --git a/opm/geomech/eclgeomechmodel.hh b/opm/geomech/eclgeomechmodel.hh index a552c8b..4756229 100644 --- a/opm/geomech/eclgeomechmodel.hh +++ b/opm/geomech/eclgeomechmodel.hh @@ -83,35 +83,32 @@ namespace Opm{ include_fracture_contributions_ = param.get("include_fracture_contributions"); //param.read("fractureparam.json"); const auto& schedule = this->simulator_.vanguard().schedule(); - int reportStepIdx = simulator_.episodeIndex(); + const int reportStepIdx = simulator_.episodeIndex(); const std::vector& wells = problem.wellModel().getLocalWells(reportStepIdx); //const std::vector& wells = schedule.getWells(reportStepIdx); //const Opm::EclipseGrid& eclgrid = simulator_.vanguard().eclState().getInputGrid(); const auto& grid = simulator_.vanguard().grid(); - std::string outputDir = Parameters::Get(); - std::string caseName = simulator_.vanguard().caseName(); + const std::string outputDir = Parameters::Get(); + const std::string caseName = simulator_.vanguard().caseName(); param.put("outputdir", outputDir); param.put("casename", caseName); - fracturemodel_ = std::make_unique(grid, - wells, - param - ); + fracturemodel_ = std::make_unique + (grid, wells, param); // not to get the reservoir properties along the well before initialising the well // most important stress - fracturemodel_->updateReservoirWellProperties(simulator_); + fracturemodel_->updateReservoirWellProperties(simulator_); // add fractures along the wells fracturemodel_->addFractures(); //fracturemodel_->updateFractureReservoirCells(grid, eclgrid); fracturemodel_->updateFractureReservoirCells(grid); - fracturemodel_->initReservoirProperties(simulator_); - fracturemodel_->updateReservoirProperties(simulator_); + fracturemodel_->initReservoirProperties(simulator_); + fracturemodel_->updateReservoirProperties(simulator_); fracturemodel_->initFractureStates(); } - // get reservoir properties on fractures - // simulator need - fracturemodel_->updateReservoirProperties(simulator_); + // get reservoir properties on fractures simulator need + fracturemodel_->updateReservoirProperties(simulator_); fracturemodel_->solve(); // copy from apply action } @@ -131,7 +128,9 @@ namespace Opm{ } - std::vector> getExtraWellIndices(std::string wellname){ + std::vector> + getExtraWellIndices(const std::string& wellname) const + { return fracturemodel_->getExtraWellIndices(wellname); } @@ -213,6 +212,7 @@ namespace Opm{ elacticitysolver_.A.getLoadVector(), elacticitysolver_.comm()); { + const auto& problem = simulator_.problem(); int num_points = simulator_.vanguard().grid().size(3); Dune::BlockVector> fixed(3*num_points); fixed = 0.0; diff --git a/opm/geomech/eclproblemgeomech.hh b/opm/geomech/eclproblemgeomech.hh index 84584ba..16870b0 100644 --- a/opm/geomech/eclproblemgeomech.hh +++ b/opm/geomech/eclproblemgeomech.hh @@ -267,7 +267,6 @@ namespace Opm{ } Parent::endStepApplyAction(); - } void addConnectionsToWell(){ @@ -294,11 +293,9 @@ namespace Opm{ //auto mapper = simulator.vanguard().cartesianMapper(); auto& wellcontainer = this->wellModel().localNonshutWells(); for (auto& wellPtr : wellcontainer) { - auto wellName = wellPtr->name(); - const auto& wellcons = geomechModel_.getExtraWellIndices(wellName); - for (const auto& cons : wellcons) { - // simple calculated with upscaling - const auto [cell, WI, depth] = cons; + const auto& wellName = wellPtr->name(); + auto& ePerfs = extra_perfs[wellName]; + for (const auto& [cell, WI, depth] : geomechModel_.getExtraWellIndices(wellName)) { // map to cartesian const auto cartesianIdx = simulator.vanguard().cartesianIndex(cell); // get ijk @@ -321,7 +318,7 @@ namespace Opm{ /*sort_value*/ -1, /*defaut sattable*/ true); connection.setCF(WI); - extra_perfs[wellName].push_back(connection); + ePerfs.push_back(connection); } } bool commit_wellstate = false;