]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Changed a few things to facilitate parallel run.
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 5 Sep 2008 17:09:19 +0000 (17:09 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 5 Sep 2008 17:09:19 +0000 (17:09 +0000)
git-svn-id: https://svn.dealii.org/trunk@16755 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-32/step-32.cc
deal.II/lac/source/trilinos_sparse_matrix.cc
deal.II/lac/source/trilinos_vector.cc

index 452a79e58a7318d5f12e189039caeb834af08477..01fbca07d02b3cd0f1124eaa9fbf99e07aebf4f2 100644 (file)
@@ -572,8 +572,8 @@ void BoussinesqFlowProblem<dim>::setup_dofs ()
   {
     stokes_dof_handler.distribute_dofs (stokes_fe);
     DoFRenumbering::Cuthill_McKee (stokes_dof_handler);
-    DoFRenumbering::component_wise (stokes_dof_handler, stokes_block_component);
     DoFRenumbering::subdomain_wise (stokes_dof_handler);
+    DoFRenumbering::component_wise (stokes_dof_handler, stokes_block_component);
     
     stokes_constraints.clear ();
     DoFTools::make_hanging_node_constraints (stokes_dof_handler,
@@ -605,16 +605,16 @@ void BoussinesqFlowProblem<dim>::setup_dofs ()
                     n_T = temperature_dof_handler.n_dofs();
 
   pcout << "Number of active cells: "
-       << triangulation.n_active_cells()
-       << " (on "
-       << triangulation.n_levels()
-       << " levels)"
-       << std::endl
-       << "Number of degrees of freedom: "
-       << n_u + n_p + n_T
-       << " (" << n_u << '+' << n_p << '+'<< n_T <<')'
-       << std::endl
-       << std::endl;
+       << triangulation.n_active_cells()
+       << " (on "
+       << triangulation.n_levels()
+       << " levels)"
+       << std::endl
+       << "Number of degrees of freedom: "
+       << n_u + n_p + n_T
+       << " (" << n_u << '+' << n_p << '+'<< n_T <<')'
+       << std::endl
+       << std::endl;
 
   stokes_partitioner.clear();
   {
@@ -628,7 +628,7 @@ void BoussinesqFlowProblem<dim>::setup_dofs ()
 
     BlockSparsityPattern stokes_sparsity_pattern (2,2);
 
-    if (trilinos_communicator.MyPID() == 0)
+    //if (trilinos_communicator.MyPID() == 0)
       {
        BlockCompressedSetSparsityPattern csp (2,2);
 
@@ -664,7 +664,7 @@ void BoussinesqFlowProblem<dim>::setup_dofs ()
 
     BlockSparsityPattern stokes_preconditioner_sparsity_pattern (2,2);
     
-    if (trilinos_communicator.MyPID() == 0)
+    //if (trilinos_communicator.MyPID() == 0)
       {
        BlockCompressedSetSparsityPattern csp (2,2);
 
@@ -693,7 +693,6 @@ void BoussinesqFlowProblem<dim>::setup_dofs ()
                                         stokes_preconditioner_sparsity_pattern);
 
   }
-  std::cout << " bla " << std::flush;
 
   temperature_partitioner = Epetra_Map (n_T, 0, trilinos_communicator);
   {
@@ -703,7 +702,7 @@ void BoussinesqFlowProblem<dim>::setup_dofs ()
 
     SparsityPattern temperature_sparsity_pattern;
 
-    if (trilinos_communicator.MyPID() == 0)
+    //if (trilinos_communicator.MyPID() == 0)
       {
        CompressedSetSparsityPattern csp (n_T, n_T);
        DoFTools::make_sparsity_pattern (temperature_dof_handler, csp);
@@ -836,6 +835,8 @@ void BoussinesqFlowProblem<dim>::assemble_stokes_system ()
 
   stokes_rhs=0;
 
+  Vector<double> old_temperatures (old_temperature_solution);
+
   QGauss<dim>   quadrature_formula(stokes_degree+2);
   QGauss<dim-1> face_quadrature_formula(stokes_degree+2);
 
@@ -879,7 +880,7 @@ void BoussinesqFlowProblem<dim>::assemble_stokes_system ()
   std::vector<SymmetricTensor<2,dim> > grads_phi_u (dofs_per_cell);
   std::vector<double>                  div_phi_u   (dofs_per_cell);
   std::vector<double>                  phi_p       (dofs_per_cell);
-
+  
   const FEValuesExtractors::Vector velocities (0);
   const FEValuesExtractors::Scalar pressure (dim);
 
@@ -898,7 +899,7 @@ void BoussinesqFlowProblem<dim>::assemble_stokes_system ()
        local_matrix = 0;
        local_rhs = 0;
   
-       temperature_fe_values.get_function_values (old_temperature_solution
+       temperature_fe_values.get_function_values (old_temperatures
                                                   old_temperature_values);
   
        for (unsigned int q=0; q<n_q_points; ++q)
@@ -1135,9 +1136,9 @@ void BoussinesqFlowProblem<dim>::assemble_temperature_system ()
        stokes_fe_values.reinit (stokes_cell);
   
        temperature_fe_values.get_function_values (old_temperature_solution,
-                                                 old_temperature_values);
+                                                  old_temperature_values);
        temperature_fe_values.get_function_values (old_old_temperature_solution,
-                                                 old_old_temperature_values);
+                                                  old_old_temperature_values);
   
        temperature_fe_values.get_function_gradients (old_temperature_solution,
                                                      old_temperature_grads);
@@ -1145,9 +1146,9 @@ void BoussinesqFlowProblem<dim>::assemble_temperature_system ()
                                                      old_old_temperature_grads);
        
        temperature_fe_values.get_function_hessians (old_temperature_solution,
-                                                   old_temperature_hessians);
+                                                    old_temperature_hessians);
        temperature_fe_values.get_function_hessians (old_old_temperature_solution,
-                                                   old_old_temperature_hessians);
+                                                    old_old_temperature_hessians);
        
        temperature_right_hand_side.value_list (temperature_fe_values.get_quadrature_points(),
                                                gamma_values);
@@ -1157,17 +1158,17 @@ void BoussinesqFlowProblem<dim>::assemble_temperature_system ()
        
        const double nu
          = compute_viscosity (old_temperature_values,
-                             old_old_temperature_values,
-                             old_temperature_grads,
-                             old_old_temperature_grads,
-                             old_temperature_hessians,
-                             old_old_temperature_hessians,
-                             present_stokes_values,
-                             gamma_values,
-                             global_u_infty,
-                             global_T_range.second - global_T_range.first,
-                             global_Omega_diameter, cell->diameter(),
-                             old_time_step);
+                              old_old_temperature_values,
+                              old_temperature_grads,
+                              old_old_temperature_grads,
+                              old_temperature_hessians,
+                              old_old_temperature_hessians,
+                              present_stokes_values,
+                              gamma_values,
+                              global_u_infty,
+                              global_T_range.second - global_T_range.first,
+                              global_Omega_diameter, cell->diameter(),
+                              old_time_step);
        
        for (unsigned int q=0; q<n_q_points; ++q)
          {
@@ -1192,47 +1193,47 @@ void BoussinesqFlowProblem<dim>::assemble_temperature_system ()
              {
                for (unsigned int i=0; i<dofs_per_cell; ++i)
                  local_rhs(i) += ((time_step + old_time_step) / old_time_step *
-                                 old_T * phi_T[i]
-                                 -
-                                 (time_step * time_step) /
-                                 (old_time_step * (time_step + old_time_step)) *
-                                 old_old_T * phi_T[i]
-                                 -
-                                 time_step *
-                                 present_u *
-                                 ((1+time_step/old_time_step) * old_grad_T
+                                  old_T * phi_T[i]
+                                  -
+                                  (time_step * time_step) /
+                                  (old_time_step * (time_step + old_time_step)) *
+                                  old_old_T * phi_T[i]
+                                  -
+                                  time_step *
+                                  present_u *
+                                  ((1+time_step/old_time_step) * old_grad_T
                                    -
                                    time_step / old_time_step * old_old_grad_T) *
-                                 phi_T[i]
-                                 -
-                                 time_step *
-                                 nu *
-                                 ((1+time_step/old_time_step) * old_grad_T
+                                  phi_T[i]
+                                  -
+                                  time_step *
+                                  nu *
+                                  ((1+time_step/old_time_step) * old_grad_T
                                    -
                                    time_step / old_time_step * old_old_grad_T) *
-                                 grad_phi_T[i]
-                                 +
-                                 time_step *
-                                 gamma_values[q] * phi_T[i])
-                                 *
-                                 temperature_fe_values.JxW(q);
+                                  grad_phi_T[i]
+                                  +
+                                  time_step *
+                                  gamma_values[q] * phi_T[i])
+                                 * 
+                                 temperature_fe_values.JxW(q);
              }
            else
              {
                for (unsigned int i=0; i<dofs_per_cell; ++i)
                  local_rhs(i) += (old_T * phi_T[i]
-                                 -
-                                 time_step *
-                                 present_u * old_grad_T * phi_T[i]
-                                 -
-                                 time_step *
-                                 nu *
-                                 old_grad_T * grad_phi_T[i]
-                                 +
-                                 time_step *
-                                 gamma_values[q] * phi_T[i])
+                                  -
+                                  time_step *
+                                  present_u * old_grad_T * phi_T[i]
+                                  -
+                                  time_step *
+                                  nu *
+                                  old_grad_T * grad_phi_T[i]
+                                  +
+                                  time_step *
+                                  gamma_values[q] * phi_T[i])
                                  *
-                                 temperature_fe_values.JxW(q);
+                                 temperature_fe_values.JxW(q);
              }
          }
        
@@ -1272,9 +1273,9 @@ void BoussinesqFlowProblem<dim>::solve ()
     gmres.solve(stokes_matrix, stokes_solution, stokes_rhs, preconditioner);
 
     pcout << "   "
-        << solver_control.last_step()
-        << " GMRES iterations for Stokes subsystem."
-        << std::endl;
+         << solver_control.last_step()
+         << " GMRES iterations for Stokes subsystem."
+         << std::endl;
 
     stokes_constraints.distribute (stokes_solution);
   }
@@ -1304,9 +1305,9 @@ void BoussinesqFlowProblem<dim>::solve ()
     temperature_constraints.distribute (temperature_solution);
 
     pcout << "   "
-        << solver_control.last_step()
-        << " CG iterations for temperature."
-        << std::endl;
+         << solver_control.last_step()
+         << " CG iterations for temperature."
+         << std::endl;
 
     double min_temperature = temperature_solution(0),
           max_temperature = temperature_solution(0);
@@ -1319,8 +1320,8 @@ void BoussinesqFlowProblem<dim>::solve ()
       }
     
     pcout << "   Temperature range: "
-        << min_temperature << ' ' << max_temperature
-        << std::endl;
+         << min_temperature << ' ' << max_temperature
+         << std::endl;
   }
 }
 
index 14075ee65e0aa99e66085ade62785d12e70b575d..9ed483cdfb4e5954331fab1113df40b04c95f385 100755 (executable)
@@ -26,10 +26,10 @@ namespace TrilinosWrappers
     SparseMatrix::const_iterator::Accessor::
     visit_present_row ()
     {
-                                     // if we are asked to visit the
-                                     // past-the-end line, then simply
-                                     // release all our caches and go on
-                                     // with life
+                                 // if we are asked to visit the
+                                 // past-the-end line, then simply
+                                 // release all our caches and go on
+                                 // with life
       if (this->a_row == matrix->m())
        {
          colnum_cache.reset ();
@@ -38,11 +38,11 @@ namespace TrilinosWrappers
          return;
        }
       
-                                     // otherwise first flush Trilinos caches
+                                 // otherwise first flush Trilinos caches
       matrix->compress ();
 
-                                     // get a representation of the present
-                                     // row
+                                 // get a representation of the present
+                                 // row
       int ncols;
       int colnums = matrix->n();
       TrilinosScalar *values = new TrilinosScalar(colnums);
@@ -52,12 +52,12 @@ namespace TrilinosWrappers
                                                  ncols, &(values[0]));
       AssertThrow (ierr == 0, ExcTrilinosError(ierr));
 
-                                     // copy it into our caches if the line
-                                     // isn't empty. if it is, then we've
-                                     // done something wrong, since we
-                                     // shouldn't have initialized an
-                                     // iterator for an empty line (what
-                                     // would it point to?)
+                                 // copy it into our caches if the
+                                 // line isn't empty. if it is, then
+                                 // we've done something wrong, since
+                                 // we shouldn't have initialized an
+                                 // iterator for an empty line (what
+                                 // would it point to?)
       Assert (ncols != 0, ExcInternalError());
       colnum_cache.reset (new std::vector<unsigned int> (colnums,
                                                        colnums+ncols));
@@ -66,17 +66,17 @@ namespace TrilinosWrappers
   }
 
 
-                                     // The constructor is actually the
-                                     // only point where we have to check
-                                     // whether we build a serial or
-                                     // a parallel Trilinos matrix.
-                                     // In the end, it does not even
-                                     // matter how many threads there
-                                     // are, but only if we use an
-                                     // MPI compiler or a standard 
-                                     // compiler. So, one thread on
-                                     // an MPI compiler will still get
-                                     // a parallel interface.
+                                 // The constructor is actually the
+                                 // only point where we have to check
+                                 // whether we build a serial or a
+                                 // parallel Trilinos matrix.
+                                 // Actually, it does not even matter
+                                 // how many threads there are, but
+                                 // only if we use an MPI compiler or
+                                 // a standard compiler. So, even one
+                                 // thread on a configuration with
+                                 // MPI will still get a parallel
+                                 // interface.
   SparseMatrix::SparseMatrix ()
                  :
 #ifdef DEAL_II_COMPILER_SUPPORTS_MPI
@@ -160,21 +160,28 @@ namespace TrilinosWrappers
   {
     unsigned int n_rows = sparsity_pattern.n_rows();
 
+                                 // TODO: As of now, assume that the
+                                 // sparsity pattern only sits at the
+                                 // zeroth processor (completely), let
+                                 // this determine the Trilinos
+                                 // sparsity pattern on that
+                                 // processor, and only broadcast the
+                                 // pattern afterwards.
     if (row_map.Comm().MyPID() == 0)
       {
        Assert (matrix->NumGlobalRows() == (int)sparsity_pattern.n_rows(),
                ExcDimensionMismatch (matrix->NumGlobalRows(),
                                      sparsity_pattern.n_rows()));
     
-                                       // Trilinos seems to have a bug for
-                                       // rectangular matrices at this point,
-                                       // so do not check for consistent 
-                                       // column numbers here.
-                                       //
-                                       //
-                                       // this bug is filed in the Sandia
-                                       // bugzilla under #4123 and should be
-                                       // fixed for version 9.0
+                                 // Trilinos seems to have a bug for
+                                 // rectangular matrices at this point,
+                                 // so do not check for consistent 
+                                 // column numbers here.
+                                 //
+                                 //
+                                 // this bug is filed in the Sandia
+                                 // bugzilla under #4123 and should be
+                                 // fixed for version 9.0
 //        Assert (matrix->NumGlobalCols() == (int)sparsity_pattern.n_cols(),
 //             ExcDimensionMismatch (matrix->NumGlobalCols(),
 //                                   sparsity_pattern.n_cols()));
@@ -197,21 +204,17 @@ namespace TrilinosWrappers
              row_indices[col] = sparsity_pattern.column_number (row, col);
            
            matrix->InsertGlobalValues(row, row_length,
-                                     &values[0], &row_indices[0]);
+                                      &values[0], &row_indices[0]);
          }
       }
+    
+    last_action = Insert;
 
-    compress();
                                  // In the end, the matrix needs to
                                  // be compressed in order to be
-                                 // really ready. However, that is
-                                 // a collective operation, so it
-                                 // has to be called on all processes
-                                 // by the user, whereas this function
-                                 // should only be used on one processor
-                                 // since our sparsity pattern data
-                                 // types are all serial as of now.
-  }
+                                 // really ready.
+    compress();
+ }
 
 
 
@@ -219,7 +222,39 @@ namespace TrilinosWrappers
   SparseMatrix::reinit (const Epetra_Map       &input_map,
                        const SparsityPattern  &sparsity_pattern)
   {
-    reinit (input_map, input_map, sparsity_pattern);
+                                 // TODO: There seems to be problem
+                                 // in Epetra when a quadratic matrix
+                                 // is initialized with both row and
+                                 // column map. Maybe find something
+                                 // more out about this...
+    //reinit (input_map, input_map, sparsity_pattern);
+    matrix.reset();
+
+    unsigned int n_rows = sparsity_pattern.n_rows();
+
+    if (row_map.Comm().MyPID() == 0)
+      {
+       Assert (input_map.NumGlobalElements() == (int)sparsity_pattern.n_rows(),
+               ExcDimensionMismatch (input_map.NumGlobalElements(),
+                                     sparsity_pattern.n_rows()));
+       Assert (input_map.NumGlobalElements() == (int)sparsity_pattern.n_cols(),
+               ExcDimensionMismatch (input_map.NumGlobalElements(),
+                                     sparsity_pattern.n_cols()));
+      }
+
+    row_map = input_map;
+    col_map = input_map;
+
+    std::vector<int> n_entries_per_row(n_rows);
+
+    for (unsigned int row=0; row<n_rows; ++row)
+      n_entries_per_row[(int)row] = sparsity_pattern.row_length(row);
+
+    matrix = std::auto_ptr<Epetra_FECrsMatrix>
+             (new Epetra_FECrsMatrix(Copy, row_map, &n_entries_per_row[0], 
+                                     false));
+
+    reinit (sparsity_pattern);
   }
 
 
@@ -234,14 +269,14 @@ namespace TrilinosWrappers
     unsigned int n_rows = sparsity_pattern.n_rows();
 
     if (row_map.Comm().MyPID() == 0)
-    {
-    Assert (input_row_map.NumGlobalElements() == (int)sparsity_pattern.n_rows(),
-           ExcDimensionMismatch (input_row_map.NumGlobalElements(),
-                                 sparsity_pattern.n_rows()));
-    Assert (input_col_map.NumGlobalElements() == (int)sparsity_pattern.n_cols(),
-           ExcDimensionMismatch (input_col_map.NumGlobalElements(),
-                                 sparsity_pattern.n_cols()));
-    }
+      {
+       Assert (input_row_map.NumGlobalElements() == (int)sparsity_pattern.n_rows(),
+               ExcDimensionMismatch (input_row_map.NumGlobalElements(),
+                                     sparsity_pattern.n_rows()));
+       Assert (input_col_map.NumGlobalElements() == (int)sparsity_pattern.n_cols(),
+               ExcDimensionMismatch (input_col_map.NumGlobalElements(),
+                                     sparsity_pattern.n_cols()));
+      }
 
     row_map = input_row_map;
     col_map = input_col_map;
@@ -253,7 +288,7 @@ namespace TrilinosWrappers
 
     matrix = std::auto_ptr<Epetra_FECrsMatrix>
              (new Epetra_FECrsMatrix(Copy, row_map, col_map,
-                                     &n_entries_per_row[0], true));
+                                     &n_entries_per_row[0], false));
 
     reinit (sparsity_pattern);
   }
@@ -267,7 +302,7 @@ namespace TrilinosWrappers
     row_map = sparse_matrix.row_map;
     col_map = sparse_matrix.col_map;
     matrix = std::auto_ptr<Epetra_FECrsMatrix>(new Epetra_FECrsMatrix(
-                             *sparse_matrix.matrix));
+                                                  *sparse_matrix.matrix));
   }
 
 
@@ -343,9 +378,9 @@ namespace TrilinosWrappers
   void
   SparseMatrix::clear ()
   {
-                                   // When we clear the matrix,
-                                   // reset the pointer and 
-                                   // generate an empty matrix.
+                                 // When we clear the matrix, reset
+                                 // the pointer and generate an
+                                 // empty matrix.
     matrix.reset();
 #ifdef DEAL_II_COMPILER_SUPPORTS_MPI
     row_map = Epetra_Map (0,0,Epetra_MpiComm(MPI_COMM_WORLD)),
@@ -362,7 +397,7 @@ namespace TrilinosWrappers
   void
   SparseMatrix::compress ()
   {
-                                   // flush buffers
+                                 // flush buffers
     int ierr;
     if (row_map.SameAs(col_map))
       ierr = matrix->GlobalAssemble ();
@@ -450,13 +485,13 @@ namespace TrilinosWrappers
        last_action = Add;
       }
 
-                                   // we have to do above actions in any
-                                   // case to be consistent with the MPI
-                                   // communication model (see the
-                                   // comments in the documentation of
-                                   // TrilinosWrappers::Vector), but we
-                                   // can save some work if the addend is
-                                   // zero
+                                 // we have to do above actions in any
+                                 // case to be consistent with the MPI
+                                 // communication model (see the
+                                 // comments in the documentation of
+                                 // TrilinosWrappers::Vector), but we
+                                 // can save some work if the addend
+                                 // is zero
     if (value == 0)
       return;
 
@@ -480,8 +515,8 @@ namespace TrilinosWrappers
     Assert (matrix->Filled()==true,
            ExcMessage("Matrix must be compressed before invoking clear_row."));
 
-                                   // Only do this on the rows
-                                   // owned locally on this processor.
+                                 // Only do this on the rows owned
+                                 // locally on this processor.
     int local_row = matrix->LRID(row);
     if (local_row >= 0)
       {
index 5405f9c9f6bdff3159a96dd6847546462c054a54..d57c9afd14366c28c80945f7ecaf90fd1785f2c7 100755 (executable)
@@ -24,8 +24,7 @@ namespace TrilinosWrappers
 {
   namespace internal
   {
-    VectorReference::operator TrilinosScalar () const // i believe useless with
-                                                      // trilinos
+    VectorReference::operator TrilinosScalar () const
     {
       Assert (index < vector.size(),
               ExcIndexRange (index, 0, vector.size()));
@@ -131,6 +130,7 @@ namespace TrilinosWrappers
 #endif
 
     vector = std::auto_ptr<Epetra_FEVector> (new Epetra_FEVector(map));
+    last_action = Insert;
   }
 
 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.