Merge pull request #376 from blattms/debug-parallel-amg-cpr

Fixes convergence issues with the parallel CPR

examples/flow_cp.cpp

@@ -263,6 +263,13 @@ try
         new_props->setSwatInitScaling(state.saturation(),pc);
     }
 
+    bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0);
+    const double *grav = use_gravity ? &gravity[0] : 0;
+
+    std::shared_ptr<DerivedGeology> geoprops;
+    std::shared_ptr<DerivedGeology> distributed_geology;
+    geoprops.reset(new Opm::DerivedGeology(*grid, *new_props, eclipseState, false, grav));
+    distributed_geology = geoprops; // copy for serial case.
     BlackoilState distributed_state;
     std::shared_ptr<Opm::BlackoilPropsAdFromDeck> distributed_props = new_props;
     Dune::CpGrid distributed_grid = *grid;
@@ -295,11 +302,16 @@ try
                                                   *distributed_props);
         grid->scatterData(state_handle);
         grid->scatterData(props_handle);
+        // Create a distributed Geology. Some values will be updated using communication
+        // below
+        distributed_geology.reset(new Opm::DerivedGeology(distributed_grid,
+                                                           *distributed_props, eclipseState,
+                                                           false, grav));
+        Opm::GeologyDataHandle geo_handle(global_grid, distributed_grid,
+                                          *geoprops, *distributed_geology);
+        grid->scatterData(geo_handle);
     }
 
-    bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0);
-    const double *grav = use_gravity ? &gravity[0] : 0;
-
     // Solver for Newton iterations.
     std::unique_ptr<NewtonIterationBlackoilInterface> fis_solver;
     
@@ -318,13 +330,11 @@ try
     // initialize variables
     simtimer.init(timeMap);
 
-    Opm::DerivedGeology geology(distributed_grid, *distributed_props, eclipseState, false, grav);
-
     std::vector<double> threshold_pressures = thresholdPressures(eclipseState, distributed_grid);
 
     SimulatorFullyImplicitBlackoil<Dune::CpGrid> simulator(param,
                                                            distributed_grid,
-                                                           geology,
+                                                           *distributed_geology,
                                                            *distributed_props,
                                                            rock_comp->isActive() ? rock_comp.get() : 0,
                                                            *fis_solver,

opm/autodiff/GeoProps.hpp

@@ -178,6 +178,8 @@ namespace Opm
         const Vector& gravityPotential() const { return gpot_   ;}
         const Vector& z()                const { return z_      ;}
         const double* gravity()          const { return gravity_;}
+        Vector&       poreVolume()             { return pvol_   ;}
+        Vector&       transmissibility()       { return trans_  ;}
 
     private:
         template <class Grid>

opm/autodiff/RedistributeDataHandles.hpp

@@ -24,6 +24,80 @@
 
 namespace Opm
 {
+
+/// \brief a data handle to distribute Derived Geology
+class GeologyDataHandle
+{
+public:
+    /// \brief type of the data we send
+    typedef double DataType;
+    /// \brief Constructor
+    /// \param sendGrid    The grid that the data is attached to when sending.
+    /// \param recvGrid    The grid that the data is attached to when receiving.
+    /// \param sendGeology The state where we will retieve the values to be sent.
+    /// \param recvGeology The state where we will store the received values.
+    GeologyDataHandle(const Dune::CpGrid& sendGrid,
+                      const Dune::CpGrid& recvGrid,
+                      const DerivedGeology& sendGeology,
+                      DerivedGeology& recvGeology)
+        : sendGrid_(sendGrid), recvGrid_(recvGrid), sendGeology_(sendGeology),
+          recvGeology_(recvGeology)
+    {}
+
+    bool fixedsize(int /*dim*/, int /*codim*/)
+    {
+        return false;
+    }
+    template<class T>
+    std::size_t size(const T& e)
+    {
+        if ( T::codimension == 0)
+        {
+            return 1 + sendGrid_.numCellFaces(e.index());
+        }
+        else
+        {
+            OPM_THROW(std::logic_error, "Data handle can only be used for elements");
+        }
+    }
+    template<class B, class T>
+    void gather(B& buffer, const T& e)
+    {
+        assert( T::codimension == 0);
+        buffer.write(sendGeology_.poreVolume()[e.index()]);
+        for ( int i=0; i< sendGrid_.numCellFaces(e.index()); ++i )
+        {
+            buffer.write(sendGeology_.transmissibility()[sendGrid_.cellFace(e.index(), i)]);
+        }
+    }
+    template<class B, class T>
+    void scatter(B& buffer, const T& e, std::size_t size)
+    {
+        assert( T::codimension == 0);
+        double val;
+        buffer.read(val);
+        recvGeology_.poreVolume()[e.index()]=val;
+        for ( int i=0; i< recvGrid_.numCellFaces(e.index()); ++i )
+        {
+            buffer.read(val);
+            recvGeology_.transmissibility()[recvGrid_.cellFace(e.index(), i)]=val;
+        }
+    }
+    bool contains(int dim, int codim)
+    {
+        return dim==3 && codim==0;
+    }
+private:
+    /// \brief The grid that the data we send is associated with.
+    const Dune::CpGrid& sendGrid_;
+    /// \brief The grid that the data we receive is associated with.
+    const Dune::CpGrid& recvGrid_;
+    /// \brief The data to send.
+    const DerivedGeology& sendGeology_;
+    /// \brief The data to receive.
+    DerivedGeology& recvGeology_;
+};
+
 /// \brief a data handle to distribute the BlackoilState
 class BlackoilStateDataHandle
 {
@@ -34,7 +108,7 @@ public:
     /// \param sendGrid   The grid that the data is attached to when sending.
     /// \param recvGrid   The grid that the data is attached to when receiving.
     /// \param sendState  The state where we will retieve the values to be sent.
-    /// \parame recvState The state where we will store the received values.
+    /// \param recvState The state where we will store the received values.
     BlackoilStateDataHandle(const Dune::CpGrid& sendGrid,
                             const Dune::CpGrid& recvGrid,
                             const BlackoilState& sendState,
@@ -73,8 +147,10 @@ public:
         buffer.write(sendState_.rv()[e.index()]);
         buffer.write(sendState_.pressure()[e.index()]);
         buffer.write(sendState_.temperature()[e.index()]);
-        buffer.write(sendState_.saturation()[e.index()]);
-
+        for ( int i=0; i<sendState_.numPhases(); ++i )
+        {
+            buffer.write(sendState_.saturation()[e.index()*sendState_.numPhases()+i]);
+        }
         for ( int i=0; i<sendGrid_.numCellFaces(e.index()); ++i )
         {
             buffer.write(sendState_.facepressure()[sendGrid_.cellFace(e.index(), i)]);
@@ -95,7 +171,7 @@ public:
         for ( int i=0; i<recvState_.numPhases(); ++i )
         {
             buffer.read(val);
-            recvState_.surfacevol()[e.index()]=val;
+            recvState_.surfacevol()[e.index()*sendState_.numPhases()+i]=val;
         }
         buffer.read(val);
         recvState_.gasoilratio()[e.index()]=val;
@@ -105,9 +181,11 @@ public:
         recvState_.pressure()[e.index()]=val;
         buffer.read(val);
         recvState_.temperature()[e.index()]=val;
-        buffer.read(val);
-        recvState_.saturation()[e.index()]=val;
-
+        for ( int i=0; i<recvState_.numPhases(); ++i )
+        {
+            buffer.read(val);
+            recvState_.saturation()[e.index()*sendState_.numPhases()+i]=val;
+        }
         for ( int i=0; i<recvGrid_.numCellFaces(e.index()); ++i )
         {
             buffer.read(val);