adding StandardWell class

copied from the old implementation, which is the starting point for the
new refactoring

CMakeLists_files.cmake

@@ -50,6 +50,7 @@ list (APPEND MAIN_SOURCE_FILES
   opm/autodiff/WellMultiSegment.cpp
   opm/autodiff/MultisegmentWells.cpp
 	opm/autodiff/WellInterface.cpp
+	opm/autodiff/StandardWell.hpp
   opm/autodiff/MissingFeatures.cpp
   opm/polymer/PolymerState.cpp
   opm/polymer/PolymerBlackoilState.cpp
@@ -241,6 +242,7 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/autodiff/StandardWells.hpp
   opm/autodiff/StandardWells_impl.hpp
 	opm/autodiff/WellInterface.hpp
+	opm/autodiff/StandardWell.cpp
   opm/autodiff/StandardWellsDense.hpp
   opm/autodiff/StandardWellsSolvent.hpp
   opm/autodiff/StandardWellsSolvent_impl.hpp

opm/autodiff/StandardWell.cpp

@@ -0,0 +1,292 @@
+/*
+  Copyright 2017 SINTEF ICT, Applied Mathematics.
+  Copyright 2017 Statoil ASA.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include "config.h"
+
+#include <opm/autodiff/StandardWell.hpp>
+
+
+
+namespace Opm
+{
+    StandardWell::
+    StandardWell(const Well* well, const size_t time_step, const Wells* wells)
+    : WellInterface(well, time_step, wells)
+    , perf_densities_(number_of_perforations_)
+    , perf_pressure_diffs_(number_of_perforations_)
+    , well_variables_(blocksize) // the number of the primary variables
+    {
+        dune_B_.setBuildMode( Mat::row_wise );
+        dune_C_.setBuildMode( Mat::row_wise );
+        inv_dune_D_.setBuildMode( Mat::row_wise );
+    }
+
+
+
+
+
+    const std::vector<double>&
+    StandardWell::
+    perfDensities() const
+    {
+        return perf_densities_;
+    }
+
+
+
+
+
+    std::vector<double>&
+    StandardWell::
+    perfDensities()
+    {
+        return perf_densities_;
+    }
+
+
+
+
+
+    const std::vector<double>&
+    StandardWell::
+    perfPressureDiffs() const
+    {
+        return perf_pressure_diffs_;
+    }
+
+
+
+
+
+    std::vector<double>&
+    StandardWell::
+    perfPressureDiffs()
+    {
+        return perf_pressure_diffs_;
+    }
+
+
+
+
+
+    void
+    StandardWell::
+    assembleWellEq(Simulator& ebos_simulator,
+                   const double dt,
+                   WellState& well_state,
+                   bool only_wells)
+    {
+    }
+
+
+
+
+
+    void StandardWell::
+    setWellVariables(const WellState& well_state)
+    {
+        const int np = number_of_phases_;
+        const int nw = well_state.bhp().size();
+        // TODO: it should be the number of primary variables
+        // TODO: this is from the old version of StandardWellsDense, it is a coincidence, 3 phases and 3 primary variables
+        // TODO: it needs to be careful.
+        // TODO: the following code has to be rewritten later for correctness purpose.
+        for (int phase = 0; phase < np; ++phase) {
+            well_variables_[phase] = 0.0;
+            well_variables_[phase].setValue(well_state.wellSolutions()[index_of_well_ + nw * phase]);
+            well_variables_[phase].setDerivative(blocksize + phase, 1.0);
+        }
+    }
+
+
+
+
+
+    StandardWell::EvalWell
+    StandardWell::
+    getBhp() const
+    {
+        const WellControls* wc = well_controls_;
+        if (well_controls_get_current_type(wc) == BHP) {
+            EvalWell bhp = 0.0;
+            const double target_rate = well_controls_get_current_target(wc);
+            bhp.setValue(target_rate);
+            return bhp;
+        } else if (well_controls_get_current_type(wc) == THP) {
+            const int control = well_controls_get_current(wc);
+            const double thp = well_controls_get_current_target(wc);
+            const double alq = well_controls_iget_alq(wc, control);
+            const int table_id = well_controls_iget_vfp(wc, control);
+            EvalWell aqua = 0.0;
+            EvalWell liquid = 0.0;
+            EvalWell vapour = 0.0;
+            EvalWell bhp = 0.0;
+            double vfp_ref_depth = 0.0;
+
+            const Opm::PhaseUsage& pu = phaseUsage();
+
+            if (active()[ Water ]) {
+                aqua = getQs(pu.phase_pos[ Water]);
+            }
+            if (active()[ Oil ]) {
+                liquid = getQs(pu.phase_pos[ Oil ]);
+            }
+            if (active()[ Gas ]) {
+                vapour = getQs(pu.phase_pos[ Gas ]);
+            }
+            if (wellType() == INJECTOR) {
+                bhp = vfp_properties_->getInj()->bhp(table_id, aqua, liquid, vapour, thp);
+                vfp_ref_depth = vfp_properties_->getInj()->getTable(table_id)->getDatumDepth();
+            } else {
+                bhp = vfp_properties_->getProd()->bhp(table_id, aqua, liquid, vapour, thp, alq);
+                vfp_ref_depth = vfp_properties_->getProd()->getTable(table_id)->getDatumDepth();
+            }
+
+            // pick the density in the top layer
+            const double rho = perf_densities_[0];
+            // TODO: not sure whether it is always correct
+            const double well_ref_depth = perf_depth_[0];
+            // const double dp = wellhelpers::computeHydrostaticCorrection(wells(), wellIdx, vfp_ref_depth, rho, gravity_);
+            const double dp = wellhelpers::computeHydrostaticCorrection(well_ref_depth, vfp_ref_depth, rho, gravity_);
+            bhp -= dp;
+            return bhp;
+        }
+
+        return well_variables_[XvarWell];
+    }
+
+
+
+
+
+    StandardWell::EvalWell
+    StandardWell::
+    getQs(const int phase) const
+    {
+        EvalWell qs = 0.0;
+
+        const WellControls* wc = well_controls_;
+        const int np = number_of_phases_;
+
+        // the target from the well controls
+        const double target = well_controls_get_current_target(wc);
+
+        // TODO: the formulation for the injectors decides it only work with single phase
+        // surface rate injection control. Improvement will be required.
+        // Injectors
+        if (wellType() == INJECTOR) {
+            // TODO: we should not rely on comp_frac anymore, it should depend on the values in distr
+            const double comp_frac = wells().comp_frac[np*wellIdx + phaseIdx];
+            if (comp_frac == 0.0) {
+                return qs;
+            }
+
+            if (well_controls_get_current_type(wc) == BHP || well_controls_get_current_type(wc) == THP) {
+                return wellVariables_[XvarWell];
+            }
+
+            // rate control
+            // TODO: if it is reservoir volume rate, it should be wrong here
+            qs.setValue(target);
+            return qs;
+        }
+
+
+        // Producers
+        if (well_controls_get_current_type(wc) == BHP || well_controls_get_current_type(wc) == THP ) {
+            return wellVariables_[XvarWell] * wellVolumeFractionScaled(phase);
+        }
+
+    }
+
+
+
+
+
+
+    StandardWell::EvalWell
+    StandardWell::
+    wellVolumeFractionScaled(const int phase) const
+    {
+        // TODO: we should be able to set the g for the well based on the control type
+        // instead of using explicit code for g all the times
+        const WellControls* wc = well_controls_;
+        if (well_controls_get_current_type(wc) == RESERVOIR_RATE) {
+            const double* distr = well_controls_get_current_distr(wc);
+            if (distr[phase] > 0.) {
+                return wellVolumeFraction(phase) / distr[phase];
+            } else {
+                // TODO: not sure why return EvalWell(0.) causing problem here
+                // Probably due to the wrong Jacobians.
+                return wellVolumeFraction(phase);
+            }
+        }
+        std::vector<double> g = {1,1,0.01};
+        return (wellVolumeFraction(phase) / g[phase]);
+    }
+
+
+
+
+
+    StandardWell::EvalWell
+    StandardWell::
+    wellVolumeFraction(const int phase) const
+    {
+        if (phase == Water) {
+            return wellVariables_[WFrac];
+        }
+
+        if (phase == Gas) {
+            return wellVariables_[GFrac];
+        }
+
+        // Oil fraction
+        EvalWell well_fraction = 1.0;
+        if (active_[Water]) {
+            well_fraction -= wellVariables_[WFrac];
+        }
+
+        if (active_[Gas]) {
+            well_fraction -= wellVariables_[GFrac];
+        }
+        return well_fraction;
+    }
+
+
+
+
+
+    StandardWell::EvalWell
+    StandardWell::
+    wellSurfaceVolumeFraction(const int phase) const
+    {
+        EvalWell sum_volume_fraction_scaled = 0.;
+        const int np = number_of_phases_;
+        for (int p = 0; p < np; ++p) {
+            sum_volume_fraction_scaled += wellVolumeFractionScaled(p);
+        }
+
+        assert(sum_volume_fraction_scaled.value() != 0.);
+
+        return wellVolumeFractionScaled(phase) / sum_volume_fraction_scaled;
+    }
+
+}

opm/autodiff/StandardWell.hpp

@@ -0,0 +1,126 @@
+/*
+  Copyright 2017 SINTEF ICT, Applied Mathematics.
+  Copyright 2017 Statoil ASA.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+
+#ifndef OPM_STANDARDWELL_HEADER_INCLUDED
+#define OPM_STANDARDWELL_HEADER_INCLUDED
+
+
+#include "config.h"
+
+#include <opm/autodiff/WellInterface.hpp>
+
+#include<dune/common/fmatrix.hh>
+#include<dune/istl/bcrsmatrix.hh>
+#include<dune/istl/matrixmatrix.hh>
+
+#include <opm/material/densead/Math.hpp>
+#include <opm/material/densead/Evaluation.hpp>
+
+namespace Opm
+{
+
+    class StandardWell: public WellInterface
+    {
+
+    public:
+        using WellInterface::Simulator;
+        using WellInterface::WellState;
+
+        // the positions of the primary variables for StandardWell
+        // there are three primary variables, the second and the third ones are F_w and F_g
+        // the first one can be total rate (G_t) or bhp, based on the control
+        enum WellVariablePositions {
+            XvarWell = 0,
+            WFrac = 1,
+            GFrac = 2
+        };
+
+        // for now, using the matrix and block version in StandardWellsDense.
+        // TODO: for bettern generality, it should contain blocksize_field and blocksize_well.
+        // They are allowed to be different and it will create four types of matrix blocks and two types of
+        // vector blocks.
+
+        /* const static int blocksize = 3;
+        typedef double Scalar;
+        typedef Dune::FieldVector<Scalar, blocksize    > VectorBlockType;
+        typedef Dune::FieldMatrix<Scalar, blocksize, blocksize > MatrixBlockType;
+        typedef Dune::BCRSMatrix <MatrixBlockType> Mat;
+        typedef Dune::BlockVector<VectorBlockType> BVector;
+        typedef DenseAd::Evaluation<double, blocksize + blocksize> EvalWell; */
+        /* using WellInterface::EvalWell;
+        using WellInterface::BVector;
+        using WellInterface::Mat;
+        using WellInterface::MatrixBlockType;
+        using WellInterface::VectorBlockType; */
+
+        StandardWell(const Well* well, const size_t time_step, const Wells* wells);
+
+        /// the densities of the fluid  in each perforation
+        virtual const std::vector<double>& perfDensities() const;
+        virtual std::vector<double>& perfDensities();
+
+        /// the pressure difference between different perforations
+        virtual const std::vector<double>& perfPressureDiffs() const;
+        virtual std::vector<double>& perfPressureDiffs();
+
+        virtual void assembleWellEq(Simulator& ebos_simulator,
+                                    const double dt,
+                                    WellState& well_state,
+                                    bool only_wells);
+
+        virtual void setWellVariables(const WellState& well_state);
+
+        EvalWell wellVolumeFractionScaled(const int phase) const;
+
+        EvalWell wellVolumeFraction(const int phase) const;
+
+        EvalWell wellSurfaceVolumeFraction(const int phase) const;
+
+    protected:
+        // densities of the fluid in each perforation
+        std::vector<double> perf_densities_;
+        // pressure drop between different perforations
+        std::vector<double> perf_pressure_diffs_;
+
+        // TODO: probably, they should be moved to the WellInterface, when
+        // we decide the template paramters.
+        // two off-diagonal matrices
+        Mat dune_B_;
+        Mat dune_C_;
+        // diagonal matrix for the well
+        Mat inv_dune_D_;
+
+        BVector res_well_;
+
+        std::vector<EvalWell> well_variables_;
+
+        // TODO: this function should be moved to the base class.
+        // while it faces chanllenges for MSWell later, since the calculation of bhp
+        // based on THP is never implemented for MSWell yet.
+        EvalWell getBhp() const;
+
+        // TODO: it is also possible to be moved to the base class.
+        EvalWell getQs(const int phase) const;
+    };
+
+}
+
+#endif // OPM_STANDARDWELL_HEADER_INCLUDED