[cig-commits] commit by buerg to /var/svn/dealii/aspect

dealii.demon at gmail.com dealii.demon at gmail.com
Fri Aug 9 09:44:09 PDT 2013


Revision 1821

Remove Stokes preconditioner matrix.

U   trunk/aspire/include/aspect/simulator.h
U   trunk/aspire/source/simulator/assembly.cc
U   trunk/aspire/source/simulator/core.cc
U   trunk/aspire/source/simulator/solver.cc


http://www.dealii.org/websvn/revision.php?repname=Aspect+Repository&path=%2F&rev=1821&peg=1821

Diff:
Modified: trunk/aspire/include/aspect/simulator.h
===================================================================
--- trunk/aspire/include/aspect/simulator.h	2013-07-30 17:17:53 UTC (rev 1820)
+++ trunk/aspire/include/aspect/simulator.h	2013-08-09 16:43:42 UTC (rev 1821)
@@ -68,7 +68,6 @@
     {
       namespace Scratch
       {
-        template <int dim>      struct StokesPreconditioner;
         template <int dim>      struct StokesSystem;
         template <int dim>      struct TemperatureSystem;
         template <int dim>      struct CompositionSystem;
@@ -76,7 +75,6 @@
 
       namespace CopyData
       {
-        template <int dim>      struct StokesPreconditioner;
         template <int dim>      struct StokesSystem;
         template <int dim>      struct TemperatureSystem;
         template <int dim>      struct CompositionSystem;
@@ -565,16 +563,6 @@
       void setup_system_matrix_composition (const IndexSet &system_partitioning);
 
       /**
-       * Set up the size and structure of the matrix used to store the
-       * elements of the matrix that is used to build the preconditioner for
-       * the system.
-       *
-       * This function is implemented in
-       * <code>source/simulator/core.cc</code>.
-       */
-      void setup_system_preconditioner (const std::vector<IndexSet> &system_partitioning);
-
-      /**
        * @}
        */
 
@@ -582,39 +570,8 @@
        * @name Functions used in the assembly of linear systems
        * @{
        */
-      /**
-       * Initiate the assembly of the preconditioner for the Stokes system.
-       *
-       * This function is implemented in
-       * <code>source/simulator/assembly.cc</code>.
-       */
-      void assemble_stokes_preconditioner ();
 
       /**
-       * Compute the integrals for the preconditioner for the Stokes system
-       * on a single cell.
-       *
-       * This function is implemented in
-       * <code>source/simulator/assembly.cc</code>.
-       */
-      void
-      local_assemble_stokes_preconditioner (const SynchronousIterators<std_cxx1x::tuple<FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>,
-                                                                                        FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>,
-                                                                                        FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> > > &cell,
-                                            internal::Assembly::Scratch::StokesPreconditioner<dim> &scratch,
-                                            internal::Assembly::CopyData::StokesPreconditioner<dim> &data);
-
-      /**
-       * Copy the contribution to the preconditioner for the Stokes system
-       * from a single cell into the global matrix that stores these elements.
-       *
-       * This function is implemented in
-       * <code>source/simulator/assembly.cc</code>.
-       */
-      void
-      copy_local_to_global_stokes_preconditioner (const internal::Assembly::CopyData::StokesPreconditioner<dim> &data);
-
-      /**
        * Compute the integrals for the Stokes matrix and right hand side
        * on a single cell.
        *
@@ -969,7 +926,6 @@
       LinearAlgebra::BlockSparseMatrix                          system_matrix_velocity;
       LinearAlgebra::SparseMatrix                               system_matrix_temperature;
       LinearAlgebra::BlockSparseMatrix                          system_matrix_composition;
-      LinearAlgebra::BlockSparseMatrix                          system_preconditioner_matrix;
 
       LinearAlgebra::BlockVector                                solution_velocity;
       LinearAlgebra::Vector                                     solution_temperature;

Modified: trunk/aspire/source/simulator/assembly.cc
===================================================================
--- trunk/aspire/source/simulator/assembly.cc	2013-07-30 17:17:53 UTC (rev 1820)
+++ trunk/aspire/source/simulator/assembly.cc	2013-08-09 16:43:42 UTC (rev 1821)
@@ -43,103 +43,8 @@
       namespace Scratch
       {
         template <int dim>
-        struct StokesPreconditioner
+        struct StokesSystem
         {
-          StokesPreconditioner (const FiniteElement<dim> &fe_velocity,
-                                const FiniteElement<dim> &fe_temperature,
-                                const FiniteElement<dim> &fe_composition,
-                                const Quadrature<dim>    &quadrature,
-                                const Mapping<dim>       &mapping,
-                                const UpdateFlags         update_flags_velocity,
-                                const UpdateFlags         update_flags_temperature,
-                                const UpdateFlags         update_flags_composition,
-                                const unsigned int        n_compositional_fields);
-          StokesPreconditioner (const StokesPreconditioner &data);
-
-          FEValues<dim>                        fe_values_velocity;
-          FEValues<dim>                        fe_values_temperature;
-          FEValues<dim>                        fe_values_composition;
-
-
-          std::vector<double>                  phi_p;
-          std::vector<Tensor<1,dim> >          phi_u;
-          std::vector<SymmetricTensor<2,dim> > symmetric_grads_phi_u;
-          std::vector<Tensor<2,dim> >          grads_phi_u;
-          std::vector<double>                  div_phi_u;
-          std::vector<Tensor<1,dim> >          velocity_values;
-
-          typename MaterialModel::Interface<dim>::MaterialModelInputs material_model_inputs;
-          typename MaterialModel::Interface<dim>::MaterialModelOutputs material_model_outputs;
-        };
-
-
-
-        template <int dim>
-        StokesPreconditioner<dim>::
-        StokesPreconditioner (const FiniteElement<dim> &fe_velocity,
-                              const FiniteElement<dim> &fe_temperature,
-                              const FiniteElement<dim> &fe_composition,
-                              const Quadrature<dim>    &quadrature,
-                              const Mapping<dim>       &mapping,
-                              const UpdateFlags         update_flags_velocity,
-                              const UpdateFlags         update_flags_temperature,
-                              const UpdateFlags         update_flags_composition,
-                              const unsigned int        n_compositional_fields)
-          :
-          fe_values_velocity (mapping, fe_velocity, quadrature,
-                              update_flags_velocity),
-          fe_values_temperature (mapping, fe_temperature, quadrature,
-                                 update_flags_temperature),
-          fe_values_composition (mapping, fe_composition, quadrature,
-                                 update_flags_composition),
-          phi_p (fe_velocity.dofs_per_cell),
-          phi_u (fe_velocity.dofs_per_cell),
-          symmetric_grads_phi_u (fe_velocity.dofs_per_cell),
-          grads_phi_u (fe_velocity.dofs_per_cell),
-          div_phi_u (fe_velocity.dofs_per_cell),
-          velocity_values (quadrature.size()),
-          material_model_inputs (quadrature.size (), n_compositional_fields),
-          material_model_outputs (quadrature.size (), n_compositional_fields)
-        {}
-
-
-
-        template <int dim>
-        StokesPreconditioner<dim>::
-        StokesPreconditioner (const StokesPreconditioner &scratch)
-          :
-          fe_values_velocity (scratch.fe_values_velocity.get_mapping (),
-                              scratch.fe_values_velocity.get_fe (),
-                              scratch.fe_values_velocity.get_quadrature (),
-                              scratch.fe_values_velocity.get_update_flags ()),
-          fe_values_temperature (scratch.fe_values_temperature.get_mapping (),
-                                 scratch.fe_values_temperature.get_fe (),
-                                 scratch.fe_values_temperature.get_quadrature (),
-                                 scratch.fe_values_temperature.get_update_flags ()),
-          fe_values_composition (scratch.fe_values_composition.get_mapping (),
-                                 scratch.fe_values_composition.get_fe (),
-                                 scratch.fe_values_composition.get_quadrature (),
-                                 scratch.fe_values_composition.get_update_flags ()),
-          phi_p (scratch.phi_p),
-          phi_u (scratch.phi_u),
-          symmetric_grads_phi_u (scratch.symmetric_grads_phi_u),
-          grads_phi_u (scratch.grads_phi_u),
-          div_phi_u (scratch.div_phi_u),
-          velocity_values (scratch.velocity_values),
-          material_model_inputs(scratch.material_model_inputs),
-          material_model_outputs(scratch.material_model_outputs)
-        {}
-
-
-
-        // We derive the StokesSystem scratch array from the
-        // StokesPreconditioner array. We do this because all the objects that
-        // are necessary for the assembly of the preconditioner are also
-        // needed for the actual matrix system and right hand side, plus some
-        // extra data that we need for the time stepping right hand side.
-        template <int dim>
-        struct StokesSystem : public StokesPreconditioner<dim>
-        {
           StokesSystem (const FiniteElement<dim> &fe_velocity,
                         const FiniteElement<dim> &fe_temperature,
                         const FiniteElement<dim> &fe_composition,
@@ -152,11 +57,21 @@
 
           StokesSystem (const StokesSystem<dim> &data);
           
+          FEValues<dim>                              fe_values_velocity;
+          FEValues<dim>                              fe_values_temperature;
+          FEValues<dim>                              fe_values_composition;
+          
+          std::vector<double>                        phi_p;
+          std::vector<Tensor<1,dim> >                phi_u;
+          std::vector<SymmetricTensor<2,dim> >       symmetric_grads_phi_u;
+          std::vector<Tensor<2,dim> >                grads_phi_u;
+          std::vector<double>                        div_phi_u;
           std::vector<double>                        temperature_values;
           std::vector<Tensor<1, dim> >               temperature_gradients;
           std::vector<std::vector<double> >          composition_values;
           std::vector<std::vector<Tensor<1, dim> > > composition_gradients;
           std::vector<Tensor<1, dim> >               pressure_gradients;
+          std::vector<Tensor<1,dim> >                velocity_values;
           std::vector<Tensor<1,dim> >                old_velocity_values;
 
           typename MaterialModel::Interface<dim>::MaterialModelInputs material_model_inputs;
@@ -177,12 +92,17 @@
                       const UpdateFlags         update_flags_composition,
                       const unsigned int        n_compositional_fields)
           :
-          StokesPreconditioner<dim> (fe_velocity, fe_temperature,
-                                     fe_composition, quadrature,
-                                     mapping, update_flags_velocity,
-                                     update_flags_temperature,
-                                     update_flags_composition,
-                                     n_compositional_fields),
+          fe_values_velocity (mapping, fe_velocity, quadrature,
+                              update_flags_velocity),
+          fe_values_temperature (mapping, fe_temperature, quadrature,
+                                 update_flags_temperature),
+          fe_values_composition (mapping, fe_composition, quadrature,
+                                 update_flags_composition),
+          phi_p (fe_velocity.dofs_per_cell),
+          phi_u (fe_velocity.dofs_per_cell),
+          symmetric_grads_phi_u (fe_velocity.dofs_per_cell),
+          grads_phi_u (fe_velocity.dofs_per_cell),
+          div_phi_u (fe_velocity.dofs_per_cell),
           temperature_values (quadrature.size ()),
           temperature_gradients (quadrature.size ()),
           composition_values (n_compositional_fields,
@@ -190,6 +110,7 @@
           composition_gradients (n_compositional_fields,
                                  std::vector<Tensor<1, dim> > (quadrature.size ())),
           pressure_gradients (quadrature.size ()),
+          velocity_values (quadrature.size()),
           old_velocity_values (quadrature.size()),
           material_model_inputs (quadrature.size (), n_compositional_fields),
           material_model_outputs (quadrature.size (), n_compositional_fields)
@@ -201,12 +122,29 @@
         StokesSystem<dim>::
         StokesSystem (const StokesSystem<dim> &scratch)
           :
-          StokesPreconditioner<dim> (scratch),
+          fe_values_velocity (scratch.fe_values_velocity.get_mapping (),
+                              scratch.fe_values_velocity.get_fe (),
+                              scratch.fe_values_velocity.get_quadrature (),
+                              scratch.fe_values_velocity.get_update_flags ()),
+          fe_values_temperature (scratch.fe_values_temperature.get_mapping (),
+                                 scratch.fe_values_temperature.get_fe (),
+                                  scratch.fe_values_temperature.get_quadrature (),
+                                scratch.fe_values_temperature.get_update_flags ()),
+          fe_values_composition (scratch.fe_values_composition.get_mapping (),
+                                 scratch.fe_values_composition.get_fe (),
+                                 scratch.fe_values_composition.get_quadrature (),
+                                 scratch.fe_values_composition.get_update_flags ()),
+          phi_p (scratch.phi_p),
+          phi_u (scratch.phi_u),
+          symmetric_grads_phi_u (scratch.symmetric_grads_phi_u),
+          grads_phi_u (scratch.grads_phi_u),
+          div_phi_u (scratch.div_phi_u),
           temperature_values (scratch.temperature_values),
           temperature_gradients (scratch.temperature_gradients),
           composition_values (scratch.composition_values),
           composition_gradients (scratch.composition_gradients),
           pressure_gradients (scratch.pressure_gradients),
+          velocity_values (scratch.velocity_values),
           old_velocity_values (scratch.old_velocity_values),
           material_model_inputs(scratch.material_model_inputs),
           material_model_outputs(scratch.material_model_outputs)
@@ -776,167 +714,11 @@
 
   template <int dim>
   void
-  Simulator<dim>::
-  local_assemble_stokes_preconditioner (const SynchronousIterators<std_cxx1x::tuple<FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>,
-                                                                                    FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>,
-                                                                                    FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> > > &cell,
-                                        internal::Assembly::Scratch::StokesPreconditioner<dim> &scratch,
-                                        internal::Assembly::CopyData::StokesPreconditioner<dim> &data)
-  {
-    const unsigned int   dofs_per_cell   = fe_velocity.dofs_per_cell;
-    const unsigned int   n_q_points      = scratch.fe_values_velocity.n_quadrature_points;
-
-    scratch.fe_values_velocity.reinit (std_cxx1x::get<0> (cell.iterators));
-    scratch.fe_values_temperature.reinit (std_cxx1x::get<1> (cell.iterators));
-    scratch.fe_values_composition.reinit (std_cxx1x::get<2> (cell.iterators));
-
-    data.local_matrix = 0;
-
-    compute_material_model_input_values (current_linearization_point_velocity,
-                                         current_linearization_point_temperature,
-                                         current_linearization_point_composition,
-                                         scratch.fe_values_velocity,
-                                         scratch.fe_values_temperature,
-                                         scratch.fe_values_composition,
-                                         scratch.material_model_inputs);
-
-    material_model->evaluate(scratch.material_model_inputs,scratch.material_model_outputs);
-
-    const bool axisymmetric_terms = parameters.use_axisymmetric_discretization;
-    
-    scratch.fe_values_velocity[introspection.extractors.velocities].get_function_values(current_linearization_point_velocity,
-                scratch.velocity_values);
-
-    for (unsigned int q=0; q<n_q_points; ++q)
-      {
-        for (unsigned int k=0; k<dofs_per_cell; ++k)
-        {
-          scratch.phi_p[k] = scratch.fe_values_velocity[introspection.extractors.pressure].value (k, q);
-          scratch.phi_u[k]       = scratch.fe_values_velocity[introspection.extractors.velocities].value (k, q);
-           scratch.grads_phi_u[k]
-             = scratch.fe_values_velocity[introspection.extractors.velocities].gradient(k,q);
-           scratch.symmetric_grads_phi_u[k]
-             = scratch.fe_values_velocity[introspection.extractors.velocities].symmetric_gradient(k,q);
-           scratch.div_phi_u[k]   = scratch.fe_values_velocity[introspection.extractors.velocities].divergence (k, q);
-        }
-
-        const double eta = scratch.material_model_outputs.viscosities[q];
-        const double density = scratch.material_model_outputs.densities[q];
-        
-        for (unsigned int i=0; i<dofs_per_cell; ++i)
-          for (unsigned int j=0; j<dofs_per_cell; ++j)
-            if (fe_velocity.system_to_component_index (i).first
-                ==
-                fe_velocity.system_to_component_index (j).first)
-              data.local_matrix(i,j) += ((2.0 * eta * (scratch.symmetric_grads_phi_u[i] * scratch.symmetric_grads_phi_u[j])
-                  + density * (scratch.phi_u[i] * (scratch.grads_phi_u[j] * scratch.velocity_values[q]))
-                  + (material_model->is_compressible()
-                     ? -2.0/3.0 * eta + parameters.stabilization_graddiv
-                         : parameters.stabilization_graddiv)
-                  * (scratch.div_phi_u[i] * scratch.div_phi_u[j])
-                  + pressure_scaling * pressure_scaling
-                  * (scratch.phi_p[i] * scratch.phi_p[j])
-                  / (eta + parameters.stabilization_graddiv)))
-               * (axisymmetric_terms ?
-                  scratch.fe_values_velocity.quadrature_point (q) (0) : 1)
-               + (axisymmetric_terms ?
-                  eta * scratch.phi_u[i][0]  * scratch.phi_u[j][0]
-                      / scratch.fe_values_velocity.quadrature_point (q) (0) : 0)
-               * scratch.fe_values_velocity.JxW (q);
-      }
-
-    std_cxx1x::get<0> (cell.iterators)->get_dof_indices (data.local_dof_indices);
-  }
-
-
-
-  template <int dim>
-  void
-  Simulator<dim>::
-  copy_local_to_global_stokes_preconditioner (const internal::Assembly::CopyData::StokesPreconditioner<dim> &data)
-  {
-    current_constraints_velocity.distribute_local_to_global (data.local_matrix,
-                                                             data.local_dof_indices,
-                                                             system_preconditioner_matrix);
-  }
-
-
-
-  template <int dim>
-  void
-  Simulator<dim>::assemble_stokes_preconditioner ()
-  {
-    system_preconditioner_matrix = 0;
-
-    const QGauss<dim> quadrature_formula(parameters.stokes_velocity_degree+1);
-
-    FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
-      cell_filter_velocity_begin (IteratorFilters::LocallyOwnedCell (), dof_handler_velocity.begin_active ());
-    FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
-      cell_filter_temperature_begin (IteratorFilters::LocallyOwnedCell (), dof_handler_temperature.begin_active ());
-    FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
-      cell_filter_composition_begin (IteratorFilters::LocallyOwnedCell (), dof_handler_composition.begin_active ());
-    typedef std_cxx1x::tuple<FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>,
-                             FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>,
-                             FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> >
-      iterators_triple;
-    SynchronousIterators<iterators_triple>
-      filtered_iterators_begin (iterators_triple (cell_filter_velocity_begin,
-                                                  cell_filter_temperature_begin,
-                                                  cell_filter_composition_begin));
-    FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
-      cell_filter_velocity_end (IteratorFilters::LocallyOwnedCell (), dof_handler_velocity.end ());
-    FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
-      cell_filter_temperature_end (IteratorFilters::LocallyOwnedCell (), dof_handler_temperature.end ());
-    FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
-      cell_filter_composition_end (IteratorFilters::LocallyOwnedCell (), dof_handler_composition.end ());
-    SynchronousIterators<iterators_triple>
-      filtered_iterators_end (iterators_triple (cell_filter_velocity_end,
-                                                cell_filter_temperature_end,
-                                                cell_filter_composition_end));
-
-    WorkStream::
-    run (filtered_iterators_begin,
-         filtered_iterators_end,
-         std_cxx1x::bind (&Simulator<dim>::
-                          local_assemble_stokes_preconditioner,
-                          this,
-                          std_cxx1x::_1,
-                          std_cxx1x::_2,
-                          std_cxx1x::_3),
-         std_cxx1x::bind (&Simulator<dim>::
-                          copy_local_to_global_stokes_preconditioner,
-                          this,
-                          std_cxx1x::_1),
-         internal::Assembly::Scratch::
-         StokesPreconditioner<dim> (fe_velocity, fe_temperature, fe_composition,
-                                    quadrature_formula,
-                                    mapping,
-                                    (update_JxW_values |
-                                     update_values |
-                                     update_gradients |
-                                     update_quadrature_points),
-                                    (update_gradients | update_values),
-                                    (update_gradients | update_values),
-                                    parameters.n_compositional_fields),
-         internal::Assembly::CopyData::
-         StokesPreconditioner<dim> (fe_velocity));
-
-    system_preconditioner_matrix.compress(VectorOperation::add);
-  }
-
-
-
-  template <int dim>
-  void
   Simulator<dim>::build_stokes_preconditioner ()
   {
     computing_timer.enter_section ("   Build Stokes preconditioner");
     pcout << "   Building Stokes preconditioner ..." << std::flush;
 
-    // first assemble the raw matrices necessary for the preconditioner
-    assemble_stokes_preconditioner ();
-
     // then extract the other information necessary to build the
     // AMG preconditioners for the A and M blocks
     std::vector<std::vector<bool> > constant_modes;
@@ -944,7 +726,7 @@
                                       introspection.component_masks.velocities,
                                       constant_modes);
     
-    LinearAlgebra::PreconditionILU * prec_Mp = new LinearAlgebra::PreconditionILU();
+    LinearAlgebra::PreconditionILU * prec_Mp = new LinearAlgebra::PreconditionILU ();
     Mp_preconditioner.reset (prec_Mp);
     TrilinosWrappers::PreconditionBlockwiseDirect * prec_A = new TrilinosWrappers::PreconditionBlockwiseDirect ();
     Amg_preconditioner.reset (prec_A);
@@ -952,7 +734,7 @@
     TrilinosWrappers::PreconditionBlockwiseDirect::AdditionalData Amg_data;
     Amg_data.overlap = 4;
     
-    prec_Mp->initialize (system_preconditioner_matrix.block (1,1));
+    prec_Mp->initialize (system_matrix_velocity.block (1,1));
     prec_A->initialize (system_matrix_velocity.block (0,0), Amg_data);
     
     pcout << std::endl;
@@ -1037,8 +819,7 @@
 //          right_hand_side += (scratch.velocity_values[q] * scratch.composition_gradients[n][q]) / scratch.composition_values[n][q];
         
         right_hand_side *= pressure_scaling * scratch.fe_values_velocity.JxW (q)
-                                            * (axisymmetric_terms ?
-                                               scratch.fe_values_velocity.quadrature_point (q) (0) : 1);
+                                            * (axisymmetric_terms ? scratch.fe_values_velocity.quadrature_point (q) (0) : 1);
         
         for (unsigned int i=0; i<dofs_per_cell; ++i)
           {
@@ -1049,10 +830,10 @@
                                              ? -2.0/3.0 * eta + parameters.stabilization_graddiv
                                                  : parameters.stabilization_graddiv)
                                           * (scratch.div_phi_u[i] * scratch.div_phi_u[j])
-                                          - (pressure_scaling *
-                                             scratch.div_phi_u[i] * scratch.phi_p[j])
-                                          - (pressure_scaling *
-                                             scratch.phi_p[i] * scratch.div_phi_u[j]))
+                                          + pressure_scaling *
+                                            (pressure_scaling * scratch.phi_p[i] * scratch.phi_p[j]
+                                             - scratch.div_phi_u[i] * scratch.phi_p[j]
+                                             - scratch.phi_p[i] * scratch.div_phi_u[j]))
                                          *
                                          (axisymmetric_terms ?
                                           scratch.fe_values_velocity.quadrature_point (q) (0) : 1)
@@ -1513,13 +1294,6 @@
 namespace aspect
 {
 #define INSTANTIATE(dim) \
-  template void Simulator<dim>::local_assemble_stokes_preconditioner (const SynchronousIterators<std_cxx1x::tuple<FilteredIterator<DoFHandler<dim>::active_cell_iterator>, \
-                                                                                                                  FilteredIterator<DoFHandler<dim>::active_cell_iterator>, \
-                                                                                                                  FilteredIterator<DoFHandler<dim>::active_cell_iterator> > > &cell, \
-                                                                      internal::Assembly::Scratch::StokesPreconditioner<dim> &scratch, \
-                                                                      internal::Assembly::CopyData::StokesPreconditioner<dim> &data); \
-  template void Simulator<dim>::copy_local_to_global_stokes_preconditioner (const internal::Assembly::CopyData::StokesPreconditioner<dim> &data); \
-  template void Simulator<dim>::assemble_stokes_preconditioner (); \
   template void Simulator<dim>::build_stokes_preconditioner (); \
   template void Simulator<dim>::local_assemble_stokes_system (const SynchronousIterators<std_cxx1x::tuple<FilteredIterator<DoFHandler<dim>::active_cell_iterator>, \
                                                                                                           FilteredIterator<DoFHandler<dim>::active_cell_iterator>, \

Modified: trunk/aspire/source/simulator/core.cc
===================================================================
--- trunk/aspire/source/simulator/core.cc	2013-07-30 17:17:53 UTC (rev 1820)
+++ trunk/aspire/source/simulator/core.cc	2013-08-09 16:43:42 UTC (rev 1821)
@@ -243,7 +243,7 @@
       }
 
     pressure_scaling = material_model->reference_viscosity() / geometry_model->length_scale();
-
+    
     // make sure that we don't have to fill every column of the statistics
     // object in each time step.
     statistics.set_auto_fill_mode(true);
@@ -553,14 +553,9 @@
                                                mpi_communicator);
     Table<2, DoFTools::Coupling> coupling (dim + 1, dim + 1);
     
-    for (unsigned int i = 0; i < dim; ++i)
-    {
-      coupling[i][dim] = DoFTools::always;
-      coupling[dim][i] = DoFTools::always;
-      
-      for (unsigned int j = 0; j < dim; ++j)
+    for (unsigned int i = 0; i <= dim; ++i)
+      for (unsigned int j = 0; j <= dim; ++j)
         coupling[i][j] = DoFTools::always;
-    }
     
     DoFTools::make_sparsity_pattern (dof_handler_velocity,
                                      coupling, sp,
@@ -627,36 +622,7 @@
   }
 
 
-
   template <int dim>
-  void Simulator<dim>::
-  setup_system_preconditioner (const std::vector<IndexSet> &system_partitioning)
-  {
-    Amg_preconditioner.reset ();
-    Mp_preconditioner.reset ();
-    T_preconditioner.reset ();
-    C_preconditioner.reset ();
-
-    system_preconditioner_matrix.clear ();
-
-    TrilinosWrappers::BlockSparsityPattern sp (system_partitioning,
-                                               mpi_communicator);
-    Table<2, DoFTools::Coupling> coupling (dim + 1, dim + 1);
-    
-    coupling[dim][dim] = DoFTools::always;
-    DoFTools::make_sparsity_pattern (dof_handler_velocity,
-                                     coupling, sp,
-                                     constraints_velocity, false,
-                                     Utilities::MPI::
-                                     this_mpi_process(mpi_communicator));
-    sp.compress();
-
-    system_preconditioner_matrix.reinit (sp);
-  }
-
-
-
-  template <int dim>
   void Simulator<dim>::setup_dofs ()
   {
     computing_timer.enter_section("Setup dof systems");
@@ -774,7 +740,6 @@
     setup_system_matrix_velocity (introspection.index_sets.velocity_partitioning);
     setup_system_matrix_temperature (introspection.index_sets.temperature_partitioning);
     setup_system_matrix_composition (introspection.index_sets.composition_partitioning);
-    setup_system_preconditioner (introspection.index_sets.velocity_partitioning);
 
     system_rhs_velocity.reinit(introspection.index_sets.velocity_partitioning, mpi_communicator);
     current_linearization_point_velocity.reinit (introspection.index_sets.velocity_relevant_partitioning, mpi_communicator);
@@ -1072,7 +1037,6 @@
             = solution_velocity.block (introspection.block_indices.velocities);
           current_linearization_point_velocity.block (introspection.block_indices.pressure)
             = solution_velocity.block (introspection.block_indices.pressure);
-          pressure_scaling_old = pressure_scaling;
           assemble_temperature_system ();
           build_temperature_preconditioner (T_preconditioner);
           solve_temperature ();

Modified: trunk/aspire/source/simulator/solver.cc
===================================================================
--- trunk/aspire/source/simulator/solver.cc	2013-07-30 17:17:53 UTC (rev 1820)
+++ trunk/aspire/source/simulator/solver.cc	2013-08-09 16:43:42 UTC (rev 1821)
@@ -174,9 +174,6 @@
          * @brief Constructor
          *
          * @param S The entire Stokes matrix
-         * @param Spre The matrix whose blocks are used in the definition of
-         *     the preconditioning of the Stokes matrix, i.e. containing approximations
-         *     of the A and S blocks.
          * @param Mppreconditioner Preconditioner object for the Schur complement,
          *     typically chosen as the mass matrix.
          * @param Apreconditioner Preconditioner object for the matrix A.
@@ -184,10 +181,9 @@
          *     the matrix $A$, or only apply one preconditioner step with it.
          **/
         BlockSchurPreconditioner (const LinearAlgebra::BlockSparseMatrix  &S,
-                                  const LinearAlgebra::BlockSparseMatrix  &Spre,
                                   const PreconditionerMp *Mppreconditioner,
                                   const PreconditionerA *Apreconditioner,
-                                  const bool                                  do_solve_A);
+                                  const bool do_solve_A);
 
         /**
          * Matrix vector product with this preconditioner object.
@@ -205,9 +201,8 @@
          * References to the various matrix object this preconditioner works on.
          */
         const LinearAlgebra::BlockSparseMatrix &stokes_matrix;
-        const LinearAlgebra::BlockSparseMatrix &stokes_preconditioner_matrix;
-        const PreconditionerMp                    *mp_preconditioner;
-        const PreconditionerA                     *a_preconditioner;
+        const PreconditionerMp *mp_preconditioner;
+        const PreconditionerA *a_preconditioner;
 
         /**
          * Whether to actually invert the $	ilde A$ part of the preconditioner matrix
@@ -220,14 +215,12 @@
     template <class PreconditionerA, class PreconditionerMp>
     BlockSchurPreconditioner<PreconditionerA, PreconditionerMp>::
     BlockSchurPreconditioner (const LinearAlgebra::BlockSparseMatrix  &S,
-                              const LinearAlgebra::BlockSparseMatrix  &Spre,
-                              const PreconditionerMp                     *Mppreconditioner,
-                              const PreconditionerA                      *Apreconditioner,
-                              const bool                                  do_solve_A)
+                              const PreconditionerMp *Mppreconditioner,
+                              const PreconditionerA *Apreconditioner,
+                              const bool do_solve_A)
       :
       iters_A (0),
       stokes_matrix     (S),
-      stokes_preconditioner_matrix     (Spre),
       mp_preconditioner (Mppreconditioner),
       a_preconditioner  (Apreconditioner),
       do_solve_A        (do_solve_A)
@@ -253,10 +246,10 @@
         // convergence without
         // iterating. We simply skip
         // solving in this case.
-        if (src.block(1).l2_norm() > 1e-50 || dst.block(1).l2_norm() > 1e-50)
-          solver.solve(stokes_preconditioner_matrix.block(1,1),
-                       dst.block(1), src.block(1),
-                       *mp_preconditioner);
+        if ((src.block (1).l2_norm () > 1e-50) || (dst.block (1).l2_norm () > 1e-50))
+          solver.solve (stokes_matrix.block (1, 1),
+                        dst.block (1), src.block (1),
+                        *mp_preconditioner);
 
         dst.block(1) *= -1.0;
       }
@@ -415,10 +408,9 @@
         // otherwise give it a try with a preconditioner that consists
         // of only a single V-cycle
         const internal::BlockSchurPreconditioner<LinearAlgebra::PreconditionBase,
-              LinearAlgebra::PreconditionBase>
-              preconditioner (system_matrix_velocity, system_preconditioner_matrix,
-                              &(*Mp_preconditioner), &(*Amg_preconditioner),
-                              false);
+                                                 LinearAlgebra::PreconditionBase>
+          preconditioner (system_matrix_velocity, &(*Mp_preconditioner), &(*Amg_preconditioner),
+                          false);
 
         SolverFGMRES<LinearAlgebra::BlockVector>
         solver(solver_control_cheap, mem,
@@ -433,10 +425,9 @@
     catch (SolverControl::NoConvergence)
       {
         const internal::BlockSchurPreconditioner<LinearAlgebra::PreconditionBase,
-              LinearAlgebra::PreconditionBase>
-              preconditioner (system_matrix_velocity, system_preconditioner_matrix,
-                              &(*Mp_preconditioner), &(*Amg_preconditioner),
-                              true);
+                                                 LinearAlgebra::PreconditionBase>
+          preconditioner (system_matrix_velocity, &(*Mp_preconditioner), &(*Amg_preconditioner),
+                          true);
 
         SolverFGMRES<LinearAlgebra::BlockVector>
         solver(solver_control_expensive, mem,


More information about the CIG-COMMITS mailing list