This should prevent misunderstandings about what the
well_index_on_proc is. It is not the well_index according to
the eclipse state (on open wells count) but the index of the
wells that are stored on this process' domain.
In the parallel run there are cases where wells perforate cells
that are neighbors of overlap/halo cells. On other process only
parts of the well are seen as perforations. These wells should be
ignored there. While the well was indeed ignored, the perforations
found where mistakenly added to the well found due not clearing the
wellperf_data[well_index]. This commit now does this clearing and
results in the right handling of wells for e.g. SPE9.
This is needed in parallel runs where the rock properties will not
be read from the deck but be communicated from a master process. Nevertheless
we need to be able to initialize the data structures with the correct
container size. In addition we need to be able to change the container values
from opm-autodiff's BlackoilPropsDataHandle.
This PR adds allow_cf to the wells structure that determine whether
crossflow is allowed or not. An extra argument is added to addWell(..)
to specify the allow_cf flag.
While hopefully not a bug it raises an exception with gcc's
libc debugging mode. Therefore we resort to using C++11's
std::vector::data instead.
The exception was rosen when running SPE9 in parallel.
This commit introduces a new public method, activeRegions(), that
retrieves those region IDs that contain at least one active cell.
We furthermore extend the cells() method to support lookup of
arbitrary region IDs. Non-active region IDs produce empty cell
ranges.
Intended use case is
for (const auto& reg : rmap.activeRegions()) {
const auto& c = rmap.cells(reg);
// use c
}
calculations
The dz calculated in WellDetails::getCubeDim is not correct in cases
where the face centroid of the horizontal faces is located above or
below the face centroid or the vertical faces. The cell thickness in
EclipseGrid, calculated using the Z-coordinates, is therefore used
instead.
Unused methods and arguments have been removed,
and we avoid including the EclMaterialLawManager
in the header. Clients IncompPropertiesFromDeck
and BlackoilPropertiesFromDeck have been updated.
If on one process a well completion is next to border then
it might also be stored in the neighbor process. Still not
all the completions of the well are known to the neighbor.
This breaks the previous assumption that for each well all
completions must belong to the partition of the process.
Therefore with this commit we allow wells that only have a
part of their completions assigned to the partition of the process.
This wells are deactivated under the assumption that they must
exist completely on another process due to the partitioning.
Previously well with just some shut completions errorneously triggered an
exception in parallel runs. This is fixed with this commit.
Due to the logic shut completions will always be marked as existing
on a process. (Initially all completions are marked as found. For
each open completion we check whether the cartesian index belongs to
the local grid. If that is not the case we mark it as not found).
Therefore we now check whether the found number of completions
is either the number of shut completions or the number of all completions.
In the former case the well is not stored on this process, and in the latter
case it is. In other cases we throw an exception.
gcc-4.7.2 (Debian 4.7.2-5) complained about:
"‘for’ loop initial declarations are only allowed in C99 mode
note: use option -std=c99 or -std=gnu99 to compile your code"
when seeing a loop like
for(int i=0; i<end; ++i)
This is fixed by moving the declaration before the for loop with
this commit. Altenatively, we could use the above option.
this makes it possible to switch to different saturation functions
again. So far the only supported function besides the default one is
the one which implements the "Stone 2" model.
this means the following changes:
- the "SatFuncGwseg" class is converted
- for now, Gwseg is the only saturation function supported by
SaturationPropsFromDeck. (will be changed in later commits.)
- the funcForCell() method of SaturationPropsFromDeck is removed as it
just occludes things
in any reasonable simulator which reads an ECL deck the deck is going
to decide which saturation function is to be used and not the outside
code. also, the table this which function will be using is not really the
calling code's business. (for any reasonable deck it is always going to
be a non-uniform table so it makes a lot of sense to avoid unnecessary
complexity IMO.)
this patch temporarily removes the ability to use anything except the
ECL default saturation function ("Gwseg"). this ability will be
restored later in this patch series.
namely BlackoilStateToFluidState which takes a BlackoilState object
and exposes it as a opm-material like fluid state object. Similar for
ExplicitArraysFluidState, which takes raw arrays.
since fluid states are a local API, the index of the cell to be used
for these two classes must be set externally. The advantage of this
concept is that it is possible to make "saturation functions" which
not only depend on saturations but also on arbitrary other quanties
(like temperature or phase composition) without having to change the
API of the "saturation" functions.
Bård spotet a bug after PR #805 was merged. Indead returning
-numeric_limits<type>::min() does not make sense for integral
values. This commit resorts to returning numeric_limits<type>::min().
Kudos to Bård for his attention.
This behaviour does not work for computing a global inner product.
Therfore this commit introduces a new function to the functor that
returns an appropriate initial value.
Previously we hardcoded float. Now we use the result_type of
the binary_function without any qualifiers. With any cv or reference
qualifiers std::numeric_limits uses a default implementation which
produces nonesense (e.g. numeric_limits<const int>::max() returns 0).
Previously, we used the setStatus method to set wells that do not
exist on the local grid to SHUT. Or at least this is what I thought
that ```well.setStatus(timestep, SHUT)```. Unfortunately, my
assumption was wrong. This was revealed while testing a parallel run
with SPE9 that threw an expeption about "Elements must be added in
weakly increasing order" in Opm::DynamicState::add(int, T). Seems like
the method name is a bit misleading.
As it turns out the WellManager has its own complete list of active
wells (shut wells are simply left out). Therefore we can use this
behaviour to our advantage: With this commit we not only exclude shut
wells from the list, but also the ones that do not exist on the local
grid. We even get rid of an ugly const_cast.
Currently, I have running a parallel SPE9 test that has not yet
aborted.
In this case the parallel index set might represent N entries (this might be the number of
cells of grid). Nevertheless, there several (n) equations/unknowns attached to each index.
In this case we construct a larger index set representing N*n unknows, where each unknown
is attached to an index.
This change only affects parallel runs.
In a parallel run each process only knows a part of the grid. Nevertheless
it does hold the complete well information. To resolve this the WellsManager
must be able to handle this case.
With this commit its constructor gets a flag indicating whether this is
a parallel run. If it is, then it does not throw if a well has cells that
are not present on the local part of the grid. Nevertheless it will check
that either all or none of the cells of a well are stored in the local part
of the grid.
Wells with no perforated cells on the local will still be present but set to SHUT.
This commit adds a verbose flag to the constructor of
SimulatorReport to allow for deactivating any
output to std:cout. This is handy for parallel runs where we only
want to print statistics on one process.
It does not make sense to report transport and pressure separately
for fully implicit solvers. It still makes sense to separate solver
from init and output though.
This commit updates the source code comment about using operator[] to
initialize the unordered map. Thanks to Bard's persistence we found
out that the cause is not the construction of the key value of type
std::string from const char* but the mapped type being a (mutable)
char* (due to C?).
This completes the PR #784.
g++-4.4 has problems converting const char* to char*
which it thinks is needed for constructing std::string.
Using operator[] circumvents this problem.
The compiler error fixed here was:
/usr/include/c++/4.4/bits/stl_pair.h: In constructor ‘std::pair<_T1, _T2>::pair(std::pair<_U1, _U2>&&) [with _U1 = const char*, _U2 = const char*, _T1 = const std::basic_string<char, std::char_traits<char>, std::allocator<char> >, _T2 = char*]’:
/home/mblatt/src/dune/opm/opm-core/opm/core/linalg/LinearSolverPetsc.cpp:40: instantiated from here
/usr/include/c++/4.4/bits/stl_pair.h:107: error: invalid conversion from ‘const char*’ to ‘char*’
make[2]: *** [CMakeFiles/opmcore.dir/opm/core/linalg/LinearSolverPetsc.cpp.o] Fehler 1
Currently the keyword EQUIL is not supported by the fully
implicit blackoil simulator when using CpGrid. This
commit is a first step towards this as it makes the
implementation of initStateEquil generic.
Well, you never know. There are containers that use a signed integer
for storing its size. This results in a warning about comparing signed with
unsigned integers. This commit prevents this by explicitly casting the size
to std::size_t.
Due to the size of the overlap layer and the discretization scheme
the rhs might not contain correct values for overlap cells. This
commit makes sure they are correct by an additional communication step.
One would think that such an assumption is safe in any case,
wouldn't one? But foen Eigen'S container this does not hold.
They do not provide STL compliant iterators, and access to them.
With this patch make the even stricter assumption that the containers
are random access and use operator[] instead of iterators.
As there are no functors for computing the minimum and maximum,
we convert the std::max and std::min function pointers to
functors (which is not really nice.) Previously we were somehow
tricked into using std::greater and std::less, which of course do
return true or false and not what we need. Additionally, do more
excessive testing with different ranges.
We need to compute quite a few global reductions in the
Newton method of opm-autodiff. This commit adds the functionality
to compute several reductions combined using only one global
communication. Compiles and test succeeds with one or more process.
With this commit the WellsManager will check the status of completions
before adding them to the internal struct wells
datastructure. Completions can be in the four states:
OPEN, SHUT, AUTO, POPN
Completions with state == SHUT will be ignored, wheras the wellsmanager
will throw if the states AUTO or POPN are encountered. The WELOPEN
keyword can also have the value 'STOP'; for completions that is
translated to 'SHUT' by Schedule object.
Rename the the meaning for shut as whats used in Eclipse.
STOP: Well stopped off above the formation. I.e. allow for flow in the
well.
SHUT: Well completely isolated from the formation. The well is removed
from the well list.
The old test was simply wrong: it computed the M-distance and compared
to the grid radius, which becomes dependent on the scaling of the
metric M. The corrected test in isClose() depends on the anisotropy
ratio of M and the grid radius.
Note that this patch does not introduce any real temperature
dependence but only changes the APIs for the viscosity and for the
density related methods. Note that I also don't like the fact that
this requires so many changes to so many files, but with the current
design of the property classes I cannot see a way to avoid this...
There were to identical if statements and the second one was followed
by an else branch. While in this case (if statement just throws) it is not
a bug, this commit cleans up one of the if statements.
gcc warned about the following
/home/mblatt/src/dune/opm/opm-core/opm/core/wells/WellsManager.cpp: In function ‘std::array<long unsigned int, 3ul> WellsManagerDetail::directionIndices(Opm::CompletionDirection::DirectionEnum)’:
/home/mblatt/src/dune/opm/opm-core/opm/core/wells/WellsManager.cpp:191: warning: control reaches end of non-void function
To calm it I introduced a throw clause after the switch statements. Thus adding a new
enum value will raise a warning on smart compilers, hopefully.
Shut wells are not added to the well list and the well index should
therefore not be increased when well control is set. This is similar to
whats is done for shut wells in createWellsFromSpecs.
Shut wells are not added to the well list and thus not considered in the
simulator.
The shut well test in test_wellsmanager is modified to assert this
behaviour.
BUG: This change provokes an assert in the EclipeWriter as number of
wells in wellstate is different from number of wells in the schedule.
The error checking macro makes it harder to read and harder to write, so
instead we now only check for functions that can contain errors. Bounds and
range checks are handled by PETSc and not OPM.
The previous implementation set plenty of values in the initialization list and
immediately overwrote these values with values looked up from the param group.
This patch makes it look up the parameteres from the param group argument,
making the constructor simpler.
Petsc only supports initialisation through the ParameterGroup constructor.
Calling the default, non-arg constructor is a static error, and not
implementing it makes using it break compiles.
call_petsc.c was really a thin C wrapper around the call to petsc itself and
turns out was mostly unnecessary C++ emulation. This removes the file entirely
and ports its functionality into LinearSolverPetsc.cpp.
All features from the file should now be more readable as well as properly
utilising modern C++ features.
The patch uses the CHKERRXX macro from petsc to handle errors reported by
petsc, and currently does not handle this and give the control back to OPM's
error/throw system.
Commit 96cf137 introduced support for Peaceman index calculation
that honoured general completion directions (X,Y,Z). This was
accomplished through a permutation index that reordered the
permeability and geometric extent components according to a local
coordinate system along the completion.
In a complete breakdown of logic, however, the d-component extent
vector was indexed as though it were a d-by-d matrix. This commit
restores sanity to the processing.
Pointy hat: @bska.
This commit extends the feature set of the WellsManager to support
horizontal ("X" and "Y") completions and include the net-to-gross
ratio in the Peaceman index ("Completion Transmissibility Factor,
CTF") of a well completion. The NTG factor is included if present
in the input deck represented by the "eclipseState".
There are two separate, though related, parts to this commit. The
first part splits the calculation of Peaceman's "effective radius"
out to a separate utility function, effectiveRadius(), and
generalises WellsManagerDetail::computeWellIndex() to account for
arbitrary directions and NTG factors. The second part uses
GridPropertyAccess::Compressed<> to extract the NTG vector from the
input if present while providing a fall-back value of 1.0 if no such
vector is available.
Note: We may wish to make the extraction policy configurable at some
point in the future.
This commit tightens the function header of method
WellsManager::createWellsFromSpecs()
to accept a reference-to-const 'cartesian_to_compressed' map. It
used to be a complete, copy-constructed object, so this is a slight
performance enhancement as we no longer need to copy a (somewhat)
large object on every call to the method.
This commit generalises the implementation of utility function
'getCubeDim' to support arbitrary number of space dimensions. In
actual practice there's no change in features as we only really use
a compile-time constant (= 3) to specify the number of space
dimensions.
This is a demonstration of using the
GridPropertyAccess::Compressed<>
class template. We save (some) memory by not creating the zero
fall-back vector in assignPermeability(), preferring instead to use
the fall-back/default mechanism of ArrayPolicy::ExtractFromDeck<>.
While here, adjust vector<PermComponent>::reserve() capacity to
reflect actual requirements.
Clients expect column-major (Fortran) ordering of the contiguous
"permeability_" array so that's what we create despite "tensor"
being row-major.
Suggested by: [at] atgeirr
This commit switches the assignment
diagonal = max(diagonal, minval)
to using a reference in the "diagonal" expression. This guarantees
that the indexing is done once which eases maintainability. While
here, replace the hard-coded dimension stride ('3') with the current
run-time dimension. This is mostly for symmetry because the overall
code is only really supported in three space dimension.
Calling code relies on permeability tensors being stored in column
major order (row index cycling the most rapidly). Honour that
requirement. The previous assignment implied row major ordering
(column index cycling the most rapidly). This, however, is a
pedantic rather than visible change because the surrounding code
enforces symmetric tensors whence both orderings produce the same
results when the array is viewed contiguously.
For constant capillar pressure function the saturation is
determined by cell depths:
Sg_max, Sw_min
----- goc ----
Sg_min, Sw_min
----- woc ----
Sg_min, Sw_max
This commit removes several instances of EOL whitespace in function
'swatInitScaling()'. Aesthetic only. No functional changes.
Suggested by: [at] atgeirr
The constant 1.0e-8 was used as a threshold to distinguish "low"
from "high" capillary pressure values. Introduce acutual constant
"pc_low_threshold" to clarify that role.
Suggested by: [at] atgeirr
This commit renames the 'np' parameter used to allocate small arrays
for saturations and capillary pressures to 'max_np' to better
reflect its purpose.
Suggested by: [at] atgeirr
New function well_controls_clone(), implemented in terms of the
public API only, mirrors the objective of function clone_wells(),
only for well control sets. Add a basic test to demonstrate the
function too.
The 'cpty' field is for internal memory management purposes only.
No client can know of its existence, let alone inspect or directly
change the value, so it should not be used to adjudicate control set
equality. This was useful during the refactoring work to introduce
the opm-parser support, but its utility has since ceased.
Okay'ed by: [at] atgeirr and [at] joakim-hove
Current version executes reordered solve once for each tracer. The benefit
is a simpler code and the ability to use MDU with tracers. The cost is
potentially higher runtime, compared to doing a single sweep for all
tracers (and tof).
When this boolean parameter is true (the default), tracer solutions
will be normalized so that the tracer averages will sum to one in each cell.
This behaviour is the same as before, the change is that it can now be
turned off.
this basically means using Opm::EclipseState instead of the raw deck
for these keywords.
with this, property modifiers like ADD, MULT, COPY and friends are
supported for at least the PERM* keywords. If additional keywords are
required these can be added relatively easily as well.
no ctest regressions have been observed with this patch on my machine.
The queried keywords are unknown to the parser and after manual
inspection also to the Eclipse RM. There might be wrong keywords some
left, but these were to ones needed to get SPE9 started using
sim_fibo_ad...
the largest change is that all classes below opm/core/props/pvt take
the PVT region index as an argument, the higher-level ones (i.e.,
BlackoilProps*) take cell indices.
Previous fix was wrong because it called getKeyword() outside the
hasKeyword() check. Current version (like original code) uses the
short-circuiting behaviour of && to ensure this.
This was missing in commit 4c2120c and produced some build failures
that were hard to analyse. I'm not convinced that the underlying
problem is solved, but this does at least restore the build.
The class OwnerOverlapCopyCommunication is not defined unless MPI is
avilable. Therefore, we cannot reference the type unless we know
that MPI is available in the current translation unit.
Removed conflicts in
opm/core/wells/WellsManager.cpp
that were due to the change
```diff
- pd.well_index = WellsManagerDetail::computeWellIndex(radius, cubical, cell_perm, completion->getDiameter());
+ pd.well_index = WellsManagerDetail::computeWellIndex(radius, cubical, cell_perm, completion->getSkinFactor());
```
in WellsManager::createWellsFromSpecs which moved from WellsManager.cpp to WellsManager_impl.hpp file in a previous commit.
As @bska pointed in my issue #518 this definition was
accidentally removed in commit 20468d1 during a merge of
upstream changes. Therefore I readd it with this patch.
- Saturations, phase pressures, and standard initialsation of RS and RV
now agree to baseline.
- Tables of RS and RV versus vertical depth (kw RSVD RVVD) have been
hardcoded for testing (need new parser) and the calculations agree to
baseline in the gas and oil zones. In the water zone there is some
differences: Our code computes saturated RS and RV using the final
phase pressures (these are modified to be consistent with saturations
and capillary pressures) while the baseline uses unmodified phase pressures.
To support this the solveSystem methods of the LinearSolverInterface gets
an optional additional template parameter of type boost::any. It can hold any
copy constructable object. In our case it is used to pass the information about
the parallelization into the solvers of dune-istl without the compiler needing to know
their type. Inside of LinearSolverIstl::solveSystem we check whether the type stored inside of
boost::any is the new ParallelIstlInformation. If this is the case we extract the information
and use the parallel solvers if available, otherwise we solve serial/sequential.
The new ParallelIstlInformation is needed as the OwnerOverlapCopyCommunication is not copy
constructable. This is indeed a design flaw that should and will fixed upstream, but for the
time being we need ParallelIstlInformation to transfer the ParallelIndexSet and RemoteIndices
objects.
This patch refactors the calls to the dune-istl solvers.
The SeqScalarProduct, and SequentialInformation is now explicitly
constructed and later used to construct the smoothers
generically. Aditionally the linear operator (MatrixAdapter) is
constructed before calling the various solver dependent solve
routines.
While this does not change the behaviour of the code it is a
preparatory step to support parallel solvers. These parallel
solvers only differ in the type of the scalarproduct and
linear operator used from the sequential ones.
This patch refactors the calls to the dune-istl solvers.
The SeqScalarProduct, and SequentialInformation is now explicitly
constructed and later used to construct the smoothers
generically. Aditionally the linear operator (MatrixAdapter) is
constructed before calling the various solver dependent solve
routines.
While this does not change the behaviour of the code it is a
preparatory step to support parallel solvers. These parallel
solvers only differ in the type of the scalarproduct and
linear operator used from the sequential ones.
This fixes further occurences of DUNE_HAS_FASTAMG that were
missed in pull request #530.
Previously we relied on the define DUNE_HAS_FAST_AMG to detect
whether these preconditioners are available. This define is only
available in the 2.2 release with cmake support. Therfore we now
addtionally test whether we are using dune-istl 2.3 or newer.
and throw an exception if "simple" is encountered...
According to Ove, gwseg should be used, because "gwseg is the model
relevant to the norne case - i.e the model eclipse uses.
The fix for the simple model has to wait for a refac of the satfunc
complex."
Previously we relied on the define DUNE_HAS_FAST_AMG to detect
whether these preconditioners are available. This define is only
available in the 2.2 release with cmake support. Therfore we now
addtionally test whether we are using dune-isl 2.3 or newer.
Patch 31c09aed was erroneous as it was trying to assing a
SaturationPropsFromDeck<SatFuncSimpleNonuniform> to a
SaturationPropsFromDeck<SatFuncSimpleUniform> in the constructor
taking the new parser. This patch fixes this.
