Added posibility to use eclipse fluid for 1D case. Corrected bug in guess for newton solver

examples/polymer_reorder.cpp

@@ -108,6 +108,7 @@ static void outputState(const UnstructuredGrid& grid,
     Opm::writeVtkData(grid, dm, vtkfile);
 
     // Write data (not grid) in Matlab format
+    dm["faceflux"] = &state.faceflux();
     for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) {
         std::ostringstream fname;
         fname << output_dir << "/" << it->first << "-" << std::setw(3) << std::setfill('0') << step << ".dat";
@@ -396,7 +397,8 @@ main(int argc, char** argv)
     const int num_transport_substeps = param.getDefault("num_transport_substeps", 1);
 
     // If we have a "deck_filename", grid and props will be read from that.
-    bool use_deck = param.has("deck_filename");
+    bool use_deck = param.getDefault("use_deck", true);
+    use_deck = param.has("deck_filename") && use_deck;
     boost::scoped_ptr<Opm::GridManager> grid;
     boost::scoped_ptr<Opm::IncompPropertiesInterface> props;
     boost::scoped_ptr<Opm::WellsManager> wells;
@@ -459,8 +461,6 @@ main(int argc, char** argv)
         wells.reset(new Opm::WellsManager());
         // Timer init.
         simtimer.init(param);
-        // Rock compressibility.
-        rock_comp.reset(new Opm::RockCompressibility(param));
         // Gravity.
         gravity[2] = param.getDefault("gravity", 0.0);
         // Init state variables (saturation and pressure).
@@ -484,6 +484,10 @@ main(int argc, char** argv)
         }
         // Init polymer properties.
         // Setting defaults to provide a simple example case.
+	bool use_deck_fluid = param.getDefault("use_deck_fluid", false);
+	if(!use_deck_fluid){
+	    // Rock compressibility.
+	    rock_comp.reset(new Opm::RockCompressibility(param));
 	    double c_max = param.getDefault("c_max_limit", 5.0);
 	    double mix_param = param.getDefault("mix_param", 1.0);
 	    double rock_density = param.getDefault("rock_density", 1000.0);
@@ -507,6 +511,12 @@ main(int argc, char** argv)
 	    polyprop.set(c_max, mix_param, rock_density, dead_pore_vol, res_factor, c_max_ads,
 			 static_cast<Opm::PolymerProperties::AdsorptionBehaviour>(ads_index),
 			 c_vals_visc,  visc_mult_vals, c_vals_ads, ads_vals);
+	}else{
+	    std::string deck_filename = param.get<std::string>("deck_filename");
+	    Opm::EclipseGridParser deck(deck_filename);
+	    rock_comp.reset(new Opm::RockCompressibility(deck));
+	    polyprop.readFromDeck(deck);
+	}
     }
 
     // Initialize polymer inflow function.

opm/polymer/TransportModelPolymer.cpp

@@ -960,6 +960,7 @@ namespace Opm
 	    mc_[cell] = mc;
 	    return;
 	}else{
+	    //*
 	    x[0] = saturation_[cell]-res[0];
 	    if((x[0]>1) || (x[0]<0)){
 		x[0] = 0.5;
@@ -968,12 +969,18 @@ namespace Opm
 	    if(x[0]>0){
 		    x[1] =  concentration_[cell]*saturation_[cell]-res[1];
 		    x[1] = x[1]/x[0];
+		    if(x[1]> polyprops_.cMax()){
+			x[1]= polyprops_.cMax()/2.0;		
+		    }
+		    if(x[1]<0){
+			x[1]=0;
+		    }
 	    }else{
 		x[1]=0;
 	    }	    	    
 	    //x[0]=0.5;x[1]=polyprops_.cMax()/2.0;
 	    res_eq.computeResidual(x, res, mc, ff);
-	    
+	    //*/
 	}