Fixing weighted gradient because of new local/global stuff

mikaem · Oct 14, 2014 · cca3ed1 · cca3ed1
1 parent fb991d3
commit cca3ed1
Show file tree

Hide file tree

Showing 5 changed files with 41 additions and 91 deletions.
diff --git a/demo/Probe.py b/demo/Probe.py
@@ -17,11 +17,10 @@
 p = Probes(x.flatten(), W)
 x = x*0.9 
 p.add_positions(x.flatten(), W)
-print 'hei'
 for i in range(6):
     p(w0)
-print 'hei'
 
 print p.array(2, "testarray")         # dump snapshot 2
 print p.array(filename="testarray")   # dump all snapshots
 print p.dump("testarray")
+
diff --git a/fenicstools/WeightedGradient.py b/fenicstools/WeightedGradient.py
@@ -3,7 +3,7 @@
 __copyright__ = "Copyright (C) 2013 " + __author__
 __license__ = "GNU Lesser GPL version 3 or any later version"
 import os, inspect
-from dolfin import compile_extension_module, Function, FunctionSpace, assemble, TrialFunction, TestFunction, dx, Matrix
+from dolfin import info, compile_extension_module, Function, FunctionSpace, assemble, TrialFunction, TestFunction, dx, Matrix
 
 fem_folder = os.path.abspath(os.path.join(inspect.getfile(inspect.currentframe()), "../fem"))
 gradient_code = open(os.path.join(fem_folder, 'gradient_weight.cpp'), 'r').read()
@@ -60,6 +60,8 @@ def weighted_gradient_matrix(mesh, i, family='CG', degree=1, constrained_domain=
             CC.append(Cp)
         return CC
     else:
-        dP = assemble(TrialFunction(S).dx(i)*TestFunction(DG)*dx)
+        dP = assemble(TrialFunction(S).dx(i)*TestFunction(DG)*dx)        
         Cp = compiled_gradient_module.compute_weighted_gradient_matrix(G, dP, dg)
+        #info(G, True)        
+        #info(dP, True)        
         return Cp
diff --git a/fenicstools/fem/gradient_weight.cpp b/fenicstools/fem/gradient_weight.cpp
@@ -26,15 +26,15 @@ namespace dolfin
     // Compute weights
     GenericVector& dg_vector = *DG.vector();  
     dg_vector.zero();
-    for (CellIterator cell(*mesh); !cell.end(); ++cell)
+    for (CellIterator cell(*mesh, "all"); !cell.end(); ++cell)
     {
       const std::vector<dolfin::la_index>& dofs
         = dofmap_u->cell_dofs(cell->index());
 
       std::fill(ws.begin(), ws.end(), 1./cell->volume());
-      dg_vector.add(ws.data(), dofs.size(), dofs.data());      
+      dg_vector.add_local(ws.data(), dofs.size(), dofs.data());    
     }  
-    dg_vector.apply("insert");  
+    dg_vector.apply("insert"); 
   }
 
   std::size_t dof_owner(std::vector<std::vector<std::size_t> > all_ranges,
@@ -55,7 +55,7 @@ namespace dolfin
     std::vector<double> values;
     std::vector<std::vector<std::size_t> > allcolumns;
     std::vector<std::vector<double> > allvalues;
-
+        
     const std::pair<std::size_t, std::size_t> row_range = A.local_range(0);
     const std::size_t m = row_range.second - row_range.first;
     GenericVector& weight = *DG.vector();
@@ -66,7 +66,7 @@ namespace dolfin
     int dm = weight_range.second-weight_range.first;
 
     const MPI_Comm mpi_comm = DG.function_space()->mesh()->mpi_comm();
-
+        
     // Communicate local_ranges of weights
     std::vector<std::vector<std::size_t> > all_ranges;
     MPI::all_gather(mpi_comm, weight_range_vec, all_ranges);
@@ -94,7 +94,7 @@ namespace dolfin
         }
       }
     }
-    
+
     // Communicate to all which weights are needed by the process
     std::vector<std::vector<std::size_t> > dofs_needed_recv;
     MPI::all_to_all(mpi_comm, dofs_needed, dofs_needed_recv);
@@ -110,7 +110,7 @@ namespace dolfin
       std::map<std::size_t, double> send_weights;
       for (std::size_t k = 0; k < dofs.size(); k++)
       {
-        weights_to_send[p].push_back(weight[dofs[k]]);
+        weights_to_send[p].push_back(weight[dofs[k]-weight_range.first]);
       }
     }
     std::vector<std::vector<double> > weights_to_send_recv;
@@ -124,7 +124,9 @@ namespace dolfin
         continue;
 
       for (std::size_t k = 0; k < dofs_needed[p].size(); k++)
-        received_weights[dofs_needed[p][k]] = weights_to_send_recv[p][k]; 
+      {
+        received_weights[dofs_needed[p][k]] = weights_to_send_recv[p][k];         
+      }
     }
 
     for (std::size_t row = 0; row < m; row++)
@@ -142,7 +144,7 @@ namespace dolfin
         }
         else
         {
-          values[i] = weight[columns[i]];  
+          values[i] = weight[columns[i]-weight_range.first];  
         }
 //        values[i] = 1./values[i];
       }
@@ -161,7 +163,7 @@ namespace dolfin
         }
         else
         {
-          w = weight[dof];  
+          w = weight[dof-weight_range.first];  
         }        
         values[i] = values[i]*w;
 //        values[i] = values[i]*values[i];
@@ -171,14 +173,15 @@ namespace dolfin
       allvalues.push_back(values);
       allcolumns.push_back(columns);
     }
+
     for (std::size_t row = 0; row < m; row++)
     {       
       // Get global row number
       const std::size_t global_row = row + row_range.first;
 
       A.setrow(global_row, allcolumns[row], allvalues[row]);
     }
-    A.apply("insert");
+    A.apply("insert");  
   }  
 
   std::shared_ptr<GenericMatrix> MatTransposeMatMult(GenericMatrix& A, GenericMatrix& B)
@@ -208,6 +211,7 @@ namespace dolfin
     Mat CC;
     PetscErrorCode ierr = MatMatMult(Ap->mat(), Bp->mat(), MAT_INITIAL_MATRIX, PETSC_DEFAULT, &CC);
     dolfin::PETScMatrix CCC = PETScMatrix(CC);
+    CCC.apply("insert");
     return CCC.copy();  
   }
 

diff --git a/tests/test_StatisticsProbes.py b/tests/test_StatisticsProbes.py
@@ -1,7 +1,7 @@
 import nose
 
 from dolfin import FunctionSpace, UnitCubeMesh, UnitSquareMesh, interpolate, \
-    Expression, VectorFunctionSpace
+    Expression, VectorFunctionSpace, MPI, mpi_comm_world
 from fenicstools import *
 from numpy import array
 
@@ -17,22 +17,10 @@ def test_StatisticsProbes_segregated_2D():
     for i in range(5):
         probes(u0, v0)
 
-    id0, p = probes[0] 
-    if len(p) > 0:     
-        assert p.number_of_evaluations() == 5
-        assert p.value_size() == 5
-
-        mean = p.mean()
-        var = p.variance()
-        nose.tools.assert_almost_equal(p[0][0], 2.5)
-        nose.tools.assert_almost_equal(p[0][4], 0.625)
-        nose.tools.assert_almost_equal(p[1][0], 0.5)
-        nose.tools.assert_almost_equal(p[1][1], 0.25)
-        nose.tools.assert_almost_equal(mean[0], 0.5)
-        nose.tools.assert_almost_equal(mean[1], 0.25)
-        nose.tools.assert_almost_equal(var[0], 0.25)
-        nose.tools.assert_almost_equal(var[1], 0.0625)
-        nose.tools.assert_almost_equal(var[2], 0.125)
+    p = probes.array()
+    if MPI.rank(mpi_comm_world()) == 0:
+        nose.tools.assert_almost_equal(p[0,0], 2.5)
+        nose.tools.assert_almost_equal(p[0,4], 0.625)
 
 def test_StatisticsProbes_segregated_3D():
     mesh = UnitCubeMesh(4, 4, 4)
@@ -47,26 +35,10 @@ def test_StatisticsProbes_segregated_3D():
     for i in range(5):
         probes(u0, v0, w0)
 
-    id0, p = probes[0] 
-    if len(p) > 0:     
-        assert p.number_of_evaluations() == 5
-        assert p.value_size() == 9
-
-        mean = p.mean()
-        var = p.variance()
-        nose.tools.assert_almost_equal(p[0][0], 2.5)
-        nose.tools.assert_almost_equal(p[0][4], 0.3125)
-        nose.tools.assert_almost_equal(p[1][0], 0.5)
-        nose.tools.assert_almost_equal(p[1][1], 0.25)
-        nose.tools.assert_almost_equal(p[1][2], 0.25)
-        nose.tools.assert_almost_equal(mean[0], 0.5)
-        nose.tools.assert_almost_equal(mean[1], 0.25)
-        nose.tools.assert_almost_equal(mean[2], 0.25)
-        nose.tools.assert_almost_equal(var[0], 0.25)
-        nose.tools.assert_almost_equal(var[1], 0.0625)
-        nose.tools.assert_almost_equal(var[2], 0.0625)
-        nose.tools.assert_almost_equal(var[3], 0.125)
-        nose.tools.assert_almost_equal(var[4], 0.125)
+    p = probes.array()
+    if MPI.rank(mpi_comm_world()) == 0:
+        nose.tools.assert_almost_equal(p[0,0], 2.5)
+        nose.tools.assert_almost_equal(p[0,4], 0.3125)
 
 def test_StatisticsProbes_vector_2D():
     mesh = UnitSquareMesh(4, 4)
@@ -79,22 +51,10 @@ def test_StatisticsProbes_vector_2D():
     for i in range(5):
         probes(u0)
 
-    id0, p = probes[0] 
-    if len(p) > 0:     
-        assert p.number_of_evaluations() == 5
-        assert p.value_size() == 5
-
-        mean = p.mean()
-        var = p.variance()
-        nose.tools.assert_almost_equal(p[0][0], 2.5)
-        nose.tools.assert_almost_equal(p[0][4], 0.625)
-        nose.tools.assert_almost_equal(p[1][0], 0.5)
-        nose.tools.assert_almost_equal(p[1][1], 0.25)
-        nose.tools.assert_almost_equal(mean[0], 0.5)
-        nose.tools.assert_almost_equal(mean[1], 0.25)
-        nose.tools.assert_almost_equal(var[0], 0.25)
-        nose.tools.assert_almost_equal(var[1], 0.0625)
-        nose.tools.assert_almost_equal(var[2], 0.125)
+    p = probes.array()
+    if MPI.rank(mpi_comm_world()) == 0:
+        nose.tools.assert_almost_equal(p[0,0], 2.5)
+        nose.tools.assert_almost_equal(p[0,4], 0.625)
 
 def test_StatisticsProbes_vector_3D():
     mesh = UnitCubeMesh(4, 4, 4)
@@ -107,26 +67,11 @@ def test_StatisticsProbes_vector_3D():
     for i in range(5):
         probes(u0)
 
-    id0, p = probes[0] 
-    if len(p) > 0:     
-        assert p.number_of_evaluations() == 5
-        assert p.value_size() == 9
+    p = probes.array()
+    if MPI.rank(mpi_comm_world()) == 0:
 
-        mean = p.mean()
-        var = p.variance()
-        nose.tools.assert_almost_equal(p[0][0], 2.5)
-        nose.tools.assert_almost_equal(p[0][4], 0.3125)
-        nose.tools.assert_almost_equal(p[1][0], 0.5)
-        nose.tools.assert_almost_equal(p[1][1], 0.25)
-        nose.tools.assert_almost_equal(p[1][2], 0.25)
-        nose.tools.assert_almost_equal(mean[0], 0.5)
-        nose.tools.assert_almost_equal(mean[1], 0.25)
-        nose.tools.assert_almost_equal(mean[2], 0.25)
-        nose.tools.assert_almost_equal(var[0], 0.25)
-        nose.tools.assert_almost_equal(var[1], 0.0625)
-        nose.tools.assert_almost_equal(var[2], 0.0625)
-        nose.tools.assert_almost_equal(var[3], 0.125)
-        nose.tools.assert_almost_equal(var[4], 0.125)
+        nose.tools.assert_almost_equal(p[0,0], 2.5)
+        nose.tools.assert_almost_equal(p[0,4], 0.3125)
 
-if __name__ == '__main__':
-    nose.run(defaultTest=__name__)
+#if __name__ == '__main__':
+    #nose.run(defaultTest=__name__)
diff --git a/tests/test_WeightedGradient.py b/tests/test_WeightedGradient.py
@@ -30,5 +30,5 @@ def test_WeightedGradient():
     du.vector()[:] = wx[2] * u.vector()
     nose.tools.assert_almost_equal(du.vector().min(), 3.)
 
-#if __name__ == '__main__':
-    #nose.run(defaultTest=__name__)
+if __name__ == '__main__':
+    nose.run(defaultTest=__name__)