Fixing weighted gradient because of new local/global stuff

mikaem · mikaem · commit cca3ed1cc392 · 2014-10-14T07:55:54.000-07:00
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__)

Original file line number	Diff line number	Diff line change
`@@ -26,15 +26,15 @@ namespace dolfin`
`26`	`26`	`// Compute weights`
`27`	`27`	`GenericVector& dg_vector = *DG.vector();`
`28`	`28`	`dg_vector.zero();`
`29`		`- for (CellIterator cell(*mesh); !cell.end(); ++cell)`
	`29`	`+ for (CellIterator cell(*mesh, "all"); !cell.end(); ++cell)`
`30`	`30`	`{`
`31`	`31`	`const std::vector<dolfin::la_index>& dofs`
`32`	`32`	`= dofmap_u->cell_dofs(cell->index());`
`33`	`33`
`34`	`34`	`std::fill(ws.begin(), ws.end(), 1./cell->volume());`
`35`		`- dg_vector.add(ws.data(), dofs.size(), dofs.data());`
	`35`	`+ dg_vector.add_local(ws.data(), dofs.size(), dofs.data());`
`36`	`36`	`}`
`37`		`- dg_vector.apply("insert");`
	`37`	`+ dg_vector.apply("insert");`
`38`	`38`	`}`
`39`	`39`
`40`	`40`	`std::size_t dof_owner(std::vector<std::vector<std::size_t> > all_ranges,`
`@@ -55,7 +55,7 @@ namespace dolfin`
`55`	`55`	`std::vector<double> values;`
`56`	`56`	`std::vector<std::vector<std::size_t> > allcolumns;`
`57`	`57`	`std::vector<std::vector<double> > allvalues;`
`58`		`-`
	`58`	`+`
`59`	`59`	`const std::pair<std::size_t, std::size_t> row_range = A.local_range(0);`
`60`	`60`	`const std::size_t m = row_range.second - row_range.first;`
`61`	`61`	`GenericVector& weight = *DG.vector();`
`@@ -66,7 +66,7 @@ namespace dolfin`
`66`	`66`	`int dm = weight_range.second-weight_range.first;`
`67`	`67`
`68`	`68`	`const MPI_Comm mpi_comm = DG.function_space()->mesh()->mpi_comm();`
`69`		`-`
	`69`	`+`
`70`	`70`	`// Communicate local_ranges of weights`
`71`	`71`	`std::vector<std::vector<std::size_t> > all_ranges;`
`72`	`72`	`MPI::all_gather(mpi_comm, weight_range_vec, all_ranges);`
`@@ -94,7 +94,7 @@ namespace dolfin`
`94`	`94`	`}`
`95`	`95`	`}`
`96`	`96`	`}`
`97`		`-`
	`97`	`+`
`98`	`98`	`// Communicate to all which weights are needed by the process`
`99`	`99`	`std::vector<std::vector<std::size_t> > dofs_needed_recv;`
`100`	`100`	`MPI::all_to_all(mpi_comm, dofs_needed, dofs_needed_recv);`
`@@ -110,7 +110,7 @@ namespace dolfin`
`110`	`110`	`std::map<std::size_t, double> send_weights;`
`111`	`111`	`for (std::size_t k = 0; k < dofs.size(); k++)`
`112`	`112`	`{`
`113`		`- weights_to_send[p].push_back(weight[dofs[k]]);`
	`113`	`+ weights_to_send[p].push_back(weight[dofs[k]-weight_range.first]);`
`114`	`114`	`}`
`115`	`115`	`}`
`116`	`116`	`std::vector<std::vector<double> > weights_to_send_recv;`
`@@ -124,7 +124,9 @@ namespace dolfin`
`124`	`124`	`continue;`
`125`	`125`
`126`	`126`	`for (std::size_t k = 0; k < dofs_needed[p].size(); k++)`
`127`		`- received_weights[dofs_needed[p][k]] = weights_to_send_recv[p][k];`
	`127`	`+ {`
	`128`	`+ received_weights[dofs_needed[p][k]] = weights_to_send_recv[p][k];`
	`129`	`+ }`
`128`	`130`	`}`
`129`	`131`
`130`	`132`	`for (std::size_t row = 0; row < m; row++)`
`@@ -142,7 +144,7 @@ namespace dolfin`
`142`	`144`	`}`
`143`	`145`	`else`
`144`	`146`	`{`
`145`		`- values[i] = weight[columns[i]];`
	`147`	`+ values[i] = weight[columns[i]-weight_range.first];`
`146`	`148`	`}`
`147`	`149`	`// values[i] = 1./values[i];`
`148`	`150`	`}`
`@@ -161,7 +163,7 @@ namespace dolfin`
`161`	`163`	`}`
`162`	`164`	`else`
`163`	`165`	`{`
`164`		`- w = weight[dof];`
	`166`	`+ w = weight[dof-weight_range.first];`
`165`	`167`	`}`
`166`	`168`	`values[i] = values[i]*w;`
`167`	`169`	`// values[i] = values[i]*values[i];`
`@@ -171,14 +173,15 @@ namespace dolfin`
`171`	`173`	`allvalues.push_back(values);`
`172`	`174`	`allcolumns.push_back(columns);`
`173`	`175`	`}`
	`176`	`+`
`174`	`177`	`for (std::size_t row = 0; row < m; row++)`
`175`	`178`	`{`
`176`	`179`	`// Get global row number`
`177`	`180`	`const std::size_t global_row = row + row_range.first;`
`178`	`181`
`179`	`182`	`A.setrow(global_row, allcolumns[row], allvalues[row]);`
`180`	`183`	`}`
`181`		`- A.apply("insert");`
	`184`	`+ A.apply("insert");`
`182`	`185`	`}`
`183`	`186`
`184`	`187`	`std::shared_ptr<GenericMatrix> MatTransposeMatMult(GenericMatrix& A, GenericMatrix& B)`
`@@ -208,6 +211,7 @@ namespace dolfin`
`208`	`211`	`Mat CC;`
`209`	`212`	`PetscErrorCode ierr = MatMatMult(Ap->mat(), Bp->mat(), MAT_INITIAL_MATRIX, PETSC_DEFAULT, &CC);`
`210`	`213`	`dolfin::PETScMatrix CCC = PETScMatrix(CC);`
	`214`	`+ CCC.apply("insert");`
`211`	`215`	`return CCC.copy();`
`212`	`216`	`}`
`213`	`217`