Merge pull request #119 from atgeirr/fix-warnings-new

Silence multiple warnings.

examples/sim_poly2p_comp_reorder.cpp

@@ -196,36 +196,6 @@ try
     bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0);
     const double *grav = use_gravity ? &gravity[0] : 0;
 
-    // Initialising src
-    int num_cells = grid->c_grid()->number_of_cells;
-    std::vector<double> src(num_cells, 0.0);
-    if (use_deck) {
-        // Do nothing, wells will be the driving force, not source terms.
-    } else {
-        // Compute pore volumes, in order to enable specifying injection rate
-        // terms of total pore volume.
-        std::vector<double> porevol;
-        if (rock_comp->isActive()) {
-            computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state.pressure(), porevol);
-        } else {
-            computePorevolume(*grid->c_grid(), props->porosity(), porevol);
-        }
-        const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0);
-        const double default_injection = use_gravity ? 0.0 : 0.1;
-        const double flow_per_sec = param.getDefault<double>("injected_porevolumes_per_day", default_injection)
-            *tot_porevol_init/unit::day;
-        src[0] = flow_per_sec;
-        src[num_cells - 1] = -flow_per_sec;
-    }
-
-    // Boundary conditions.
-    FlowBCManager bcs;
-    if (param.getDefault("use_pside", false)) {
-        int pside = param.get<int>("pside");
-        double pside_pressure = param.get<double>("pside_pressure");
-        bcs.pressureSide(*grid->c_grid(), FlowBCManager::Side(pside), pside_pressure);
-    }
-
     // Linear solver.
     LinearSolverFactory linsolver(param);
 
@@ -262,8 +232,6 @@ try
                                                rock_comp->isActive() ? rock_comp.get() : 0,
                                                wells,
                                                polymer_inflow,
-                                               src,
-                                               bcs.c_bcs(),
                                                linsolver,
                                                grav);
         SimulatorTimer simtimer;
@@ -325,8 +293,6 @@ try
                                                    rock_comp->isActive() ? rock_comp.get() : 0,
                                                    wells,
                                                    *polymer_inflow,
-                                                   src,
-                                                   bcs.c_bcs(),
                                                    linsolver,
                                                    grav);
             if (reportStepIdx == 0) {

opm/polymer/SimulatorCompressiblePolymer.cpp

@@ -95,8 +95,6 @@ namespace Opm
              const RockCompressibility* rock_comp_props,
              WellsManager& wells_manager,
              const PolymerInflowInterface& polymer_inflow,
-             const std::vector<double>& src,
-             const FlowBoundaryConditions* bcs,
              LinearSolverInterface& linsolver,
              const double* gravity);
 
@@ -146,13 +144,11 @@ namespace Opm
                                                                const RockCompressibility* rock_comp_props,
                                                                WellsManager& wells_manager,
                                                                const PolymerInflowInterface& polymer_inflow,
-                                                               const std::vector<double>& src,
-                                                               const FlowBoundaryConditions* bcs,
                                                                LinearSolverInterface& linsolver,
                                                                const double* gravity)
     {
         pimpl_.reset(new Impl(param, grid, props, poly_props, rock_comp_props,
-                              wells_manager, polymer_inflow, src, bcs, linsolver, gravity));
+                              wells_manager, polymer_inflow, linsolver, gravity));
     }
 
 
@@ -177,8 +173,6 @@ namespace Opm
                                              const RockCompressibility* rock_comp_props,
                                              WellsManager& wells_manager,
                                              const PolymerInflowInterface& polymer_inflow,
-                                             const std::vector<double>& src,
-                                             const FlowBoundaryConditions* bcs,
                                              LinearSolverInterface& linsolver,
                                              const double* gravity)
         : grid_(grid),

opm/polymer/SimulatorCompressiblePolymer.hpp

@@ -67,8 +67,6 @@ namespace Opm
         /// \param[in] rock_comp_props  if non-null, rock compressibility properties
         /// \param[in] wells_manager    well manager, may manage no (null) wells
         /// \param[in] polymer_inflow   polymer inflow controls
-        /// \param[in] src              source terms
-        /// \param[in] bcs              boundary conditions, treat as all noflow if null
         /// \param[in] linsolver        linear solver
         /// \param[in] gravity          if non-null, gravity vector
         SimulatorCompressiblePolymer(const parameter::ParameterGroup& param,
@@ -78,8 +76,6 @@ namespace Opm
                                      const RockCompressibility* rock_comp_props,
                                      WellsManager& wells_manager,
                                      const PolymerInflowInterface& polymer_inflow,
-                                     const std::vector<double>& src,
-                                     const FlowBoundaryConditions* bcs,
                                      LinearSolverInterface& linsolver,
                                      const double* gravity);
 

opm/polymer/TransportSolverTwophasePolymer.cpp

@@ -116,13 +116,6 @@ namespace
 	return std::max(std::abs(res[0]), std::abs(res[1]));
     }
 
-    bool solveNewtonStepSC(const double* , const Opm::TransportSolverTwophasePolymer::ResidualEquation&,
-			   const double*, double*);
-    bool solveNewtonStepC(const double* , const Opm::TransportSolverTwophasePolymer::ResidualEquation&,
-                          const double*, double*);
-
-
-
     // Define a piecewise linear curve along which we will look for zero of the "s" or "r" residual.
     // The curve starts at "x", goes along the direction "direction" until it hits the boundary of the box of
     // admissible values for "s" and "x" (which is given by "[x_min[0], x_max[0]]x[x_min[1], x_max[1]]").
@@ -1554,63 +1547,6 @@ namespace
 	return res_eq_.computeResidualC(x_c);
     }
 
-    bool solveNewtonStepSC(const double* xx, const  Opm::TransportSolverTwophasePolymer::ResidualEquation& res_eq,
-                           const double* res, double* x_new) {
-
-    	double dres_s_dsdc[2];
-    	double dres_c_dsdc[2];
-	double dx=1e-8;
-	double x[2];
-	if(!(x[0]>0)){
-	    x[0] = dx;
-            x[1] = 0;
-	} else {
-	    x[0] = xx[0];
-            x[1] = xx[1];
-	}
-	res_eq.computeJacobiRes(x, dres_s_dsdc, dres_c_dsdc);
-	
-	double dFx_dx=(dres_s_dsdc[0]-x[1]*dres_s_dsdc[1]);
-	double dFx_dy= (dres_s_dsdc[1]/x[0]);
-	double dFy_dx=(dres_c_dsdc[0]-x[1]*dres_c_dsdc[1]);
-	double dFy_dy= (dres_c_dsdc[1]/x[0]);
-
-	double det = dFx_dx*dFy_dy - dFy_dx*dFx_dy;
-    	if (std::abs(det) < 1e-8) {
-    	    return false;
-    	} else {
-	    x_new[0] = x[0] - (res[0]*dFy_dy - res[1]*dFx_dy)/det;
-	    x_new[1] = x[1]*x[0] - (res[1]*dFx_dx - res[0]*dFy_dx)/det;
-	    x_new[1] = (x_new[0]>0) ?  x_new[1]/x_new[0] : 0.0;
-	    return true;
-	}
-    }
-    bool solveNewtonStepC(const double* xx, const  Opm::TransportSolverTwophasePolymer::ResidualEquation& res_eq,
-                          const double* res, double* x_new) {
-	
-	double dres_s_dsdc[2];
-	double dres_c_dsdc[2];
-        double dx=1e-8;
-        double x[2];
-        if(!(x[0]>0)){
-            x[0] = dx;
-            x[1] = 0;
-        } else {
-	    x[0] = xx[0];
-            x[1] = xx[1];
-        }
-        res_eq.computeJacobiRes(x, dres_s_dsdc, dres_c_dsdc);
-        double det = dres_s_dsdc[0]*dres_c_dsdc[1] - dres_c_dsdc[0]*dres_s_dsdc[1];
-	if (std::abs(det) < 1e-8) {
-    	    return false;
-    	} else {
-    	    x_new[0] = x[0] - (res[0]*dres_c_dsdc[1] - res[1]*dres_s_dsdc[1])/det;
-	    x_new[1] = x[1] - (res[1]*dres_s_dsdc[0] - res[0]*dres_c_dsdc[0])/det;
-	    return true;
-	}
-    }
-
-
 } // Anonymous namespace
 
 

opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp

@@ -354,7 +354,7 @@ namespace Opm {
                 const ADB mu = fluidViscosity(canonicalPhaseIdx, phasePressure, state.temperature, state.rs, state.rv, cond, cells_);
                 const ADB cmax = ADB::constant(cmax_, state.concentration.blockPattern());
                 const ADB mc = computeMc(state);
-                const ADB krw_eff = polymer_props_ad_.effectiveRelPerm(state.concentration, cmax, kr, state.saturation[actph]);
+                const ADB krw_eff = polymer_props_ad_.effectiveRelPerm(state.concentration, cmax, kr);
                 const ADB inv_wat_eff_visc = polymer_props_ad_.effectiveInvWaterVisc(state.concentration, mu.value().data());
                 // Reduce mobility of water phase by relperm reduction and effective viscosity increase.
                 rq_[actph].mob = tr_mult * krw_eff * inv_wat_eff_visc;
@@ -404,7 +404,7 @@ namespace Opm {
 
             for ( int idx=0; idx<MaxNumPhases+1; ++idx )
             {
-                if (idx == MaxNumPhases || active_[idx]) { // Dealing with polymer *or* an active phase.
+                if ((idx == MaxNumPhases && has_polymer_) || active_[idx]) { // Dealing with polymer *or* an active phase.
                     auto values     = std::tuple<double,double,double>(0.0 ,0.0 ,0.0);
                     auto containers = std::make_tuple(B.col(idx),
                                                       tempV.col(idx),
@@ -430,7 +430,6 @@ namespace Opm {
             }
             info.communicator().max(&maxNormWell[0], MaxNumPhases+1);
             // Compute pore volume
-            #warning Missing polymer code for MPI version
             return std::get<1>(nc_and_pv);
         }
         else
@@ -438,9 +437,9 @@ namespace Opm {
         {
             for ( int idx=0; idx<MaxNumPhases+1; ++idx )
             {
-                if (idx == MaxNumPhases || active_[idx]) { // Dealing with polymer *or* an active phase.
+                if ((idx == MaxNumPhases && has_polymer_) || active_[idx]) { // Dealing with polymer *or* an active phase.
                     B_avg[idx] = B.col(idx).sum()/nc;
-                    maxCoeff[idx]=tempV.col(idx).maxCoeff();
+                    maxCoeff[idx] = tempV.col(idx).maxCoeff();
                     R_sum[idx] = R.col(idx).sum();
                 }
                 else
@@ -454,11 +453,6 @@ namespace Opm {
                     }
                 }
             }
-            if (has_polymer_) {
-                B_avg[MaxNumPhases] = B.col(MaxNumPhases).sum()/nc;
-                maxCoeff[MaxNumPhases]=tempV.col(MaxNumPhases).maxCoeff();
-                R_sum[MaxNumPhases] = R.col(MaxNumPhases).sum();
-            }
             // Compute total pore volume
             return geo_.poreVolume().sum();
         }

opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp

@@ -407,8 +407,8 @@ namespace {
         // Initial saturation.
         assert (not x.saturation().empty());
         const DataBlock s = Eigen::Map<const DataBlock>(& x.saturation()[0], nc, np);
-        const V sw = s.col(0);
+        const V sw0 = s.col(0);
-        vars0.push_back(sw);
+        vars0.push_back(sw0);
 
         // Initial concentration.
         assert (not x.concentration().empty());
@@ -542,7 +542,7 @@ namespace {
         const V trans = subset(geo_.transmissibility(), ops_.internal_faces);
         const std::vector<ADB> kr = computeRelPerm(state);
 		const ADB cmax = ADB::constant(cmax_, state.concentration.blockPattern());
-        const ADB krw_eff = polymer_props_ad_.effectiveRelPerm(state.concentration, cmax, kr[0], state.saturation[0]);
+        const ADB krw_eff = polymer_props_ad_.effectiveRelPerm(state.concentration, cmax, kr[0]);
         const ADB mc = computeMc(state);
         computeMassFlux(trans, mc, kr[1], krw_eff, state);
         residual_.material_balance_eq[0] = pvdt*(rq_[0].accum[1] - rq_[0].accum[0])
@@ -803,27 +803,6 @@ namespace {
 
 
 
-    std::vector<ADB>
-    FullyImplicitCompressiblePolymerSolver::computeRelPermWells(const SolutionState& state,
-                                                     const DataBlock& well_s,
-                                                     const std::vector<int>& well_cells) const
-    {
-        const int nw = wells_.number_of_wells;
-        const int nperf = wells_.well_connpos[nw];
-        const std::vector<int>& bpat = state.pressure.blockPattern();
-
-        const ADB null = ADB::constant(V::Zero(nperf), bpat);
-
-        const ADB sw = state.saturation[0];
-        const ADB so = state.saturation[1];
-        const ADB sg = null;
-        return fluid_.relperm(sw, so, sg, well_cells);
-    }
-
-
-
-
-
     std::vector<ADB>
     FullyImplicitCompressiblePolymerSolver::
 	computePressures(const SolutionState& state) const

opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp

@@ -178,11 +178,7 @@ namespace Opm {
         computeRelPerm(const SolutionState& state) const;
 
         std::vector<ADB>
-        computeRelPermWells(const SolutionState& state,
+        computePressures(const SolutionState& state) const;
-                            const DataBlock& well_s,
-                            const std::vector<int>& well_cells) const;
-		std::vector<ADB>
-		computePressures(const SolutionState& state) const;
 
         void
         computeMassFlux(const int               actph ,

opm/polymer/fullyimplicit/PolymerPropsAd.cpp

@@ -244,8 +244,7 @@ namespace Opm {
     ADB
     PolymerPropsAd::effectiveRelPerm(const ADB& c,
                                      const ADB& cmax_cells,
-                                     const ADB& krw,
+                                     const ADB& krw) const
-                                     const ADB& sw) const
     {
         const int nc = c.value().size();
         V one = V::Ones(nc);

opm/polymer/fullyimplicit/PolymerPropsAd.hpp

@@ -101,7 +101,7 @@ namespace Opm {
 		/// \param[in] relperm			Array of n relative water relperm values.
 		/// \return						Array of n adsorption values.
         ADB
-        effectiveRelPerm(const ADB& c, const ADB& cmax_cells, const ADB& krw, const ADB& sw) const;
+        effectiveRelPerm(const ADB& c, const ADB& cmax_cells, const ADB& krw) const;
 
     private:
         const PolymerProperties& polymer_props_;

opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.hpp

@@ -93,18 +93,18 @@ namespace Opm
     public:
         /// Initialise from parameters and objects to observe.
         SimulatorFullyImplicitCompressiblePolymer(const parameter::ParameterGroup& param,
-                        		                  const GridT& grid,
+                                                  const GridT& grid,
                                                   const DerivedGeology& geo,
-                                   				  BlackoilPropsAdInterface& props,
+                                                  BlackoilPropsAdInterface& props,
-                                       			  const PolymerPropsAd&    polymer_props,
+                                                  const PolymerPropsAd&    polymer_props,
-                                       			  const RockCompressibility* rock_comp_props,
+                                                  const RockCompressibility* rock_comp_props,
                                                   std::shared_ptr<EclipseState> eclipse_state,
                                                   BlackoilOutputWriter& output_writer,
                                                   Opm::DeckConstPtr& deck,
-                                       			  NewtonIterationBlackoilInterface& linsolver,
+                                                  NewtonIterationBlackoilInterface& linsolver,
-                                       			  const double* gravity);
+                                                  const double* gravity);
 
-        Solver* createSolver(const Wells* wells);
+        std::unique_ptr<Solver> createSolver(const Wells* wells);
 
         void handleAdditionalWellInflow(SimulatorTimer& timer,
                                         WellsManager& wells_manager,

opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer_impl.hpp

@@ -58,15 +58,15 @@ SimulatorFullyImplicitCompressiblePolymer(const parameter::ParameterGroup& param
 template <class GridT>
 auto SimulatorFullyImplicitCompressiblePolymer<GridT>::
 createSolver(const Wells* wells)
-    -> Solver*
+    -> std::unique_ptr<Solver>
 {
-    return new Solver(BaseType::grid_,
+    return std::unique_ptr<Solver>(new Solver(BaseType::grid_,
-                      BaseType::props_,
+                                              BaseType::props_,
-                      BaseType::geo_,
+                                              BaseType::geo_,
-                      BaseType::rock_comp_props_,
+                                              BaseType::rock_comp_props_,
-                      polymer_props_,
+                                              polymer_props_,
-                      *wells,
+                                              *wells,
-                      BaseType::solver_);
+                                              BaseType::solver_));
 }
 
 template <class GridT>