]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Remove fragile schur_complement_05 test. 5586/head
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Tue, 5 Dec 2017 17:18:56 +0000 (18:18 +0100)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Tue, 5 Dec 2017 17:18:56 +0000 (18:18 +0100)
tests/lac/schur_complement_05.cc [deleted file]
tests/lac/schur_complement_05.with_trilinos=true.with_mpi=true.mpirun=1.output [deleted file]
tests/lac/schur_complement_05.with_trilinos=true.with_mpi=true.mpirun=3.output [deleted file]

diff --git a/tests/lac/schur_complement_05.cc b/tests/lac/schur_complement_05.cc
deleted file mode 100644 (file)
index 70e6bdf..0000000
+++ /dev/null
@@ -1,957 +0,0 @@
-/* ---------------------------------------------------------------------
- *
- * Copyright (C) 2008 - 2017 by the deal.II authors
- *
- * This file is part of the deal.II library.
- *
- * The deal.II library is free software; you can use it, redistribute
- * it, and/or modify it under the terms of the GNU Lesser General
- * Public License as published by the Free Software Foundation; either
- * version 2.1 of the License, or (at your option) any later version.
- * The full text of the license can be found in the file LICENSE at
- * the top level of the deal.II distribution.
- *
- * ---------------------------------------------------------------------
- */
-
-//
-// Check the operation of the Schur Complement for Trilinos vectors
-// and in parallel
-// This is a lightly modified version of the test "mpi/periodicity_02.cc"
-// The only difference is the linear solver routine, which now uses the
-// schur complement operation first tested in "lac/schur_complement_03.cc"
-//
-
-#define PERIODIC
-
-#include "../tests.h"
-#include <deal.II/base/conditional_ostream.h>
-
-#include <deal.II/distributed/grid_refinement.h>
-
-#include <deal.II/lac/solver_cg.h>
-#include <deal.II/lac/constraint_matrix.h>
-
-#include <deal.II/lac/block_sparsity_pattern.h>
-#include <deal.II/lac/trilinos_solver.h>
-#include <deal.II/lac/trilinos_precondition.h>
-#include <deal.II/lac/trilinos_block_sparse_matrix.h>
-#include <deal.II/lac/trilinos_parallel_block_vector.h>
-#include <deal.II/lac/linear_operator.h>
-#include <deal.II/lac/trilinos_linear_operator.h>
-#include <deal.II/lac/packaged_operation.h>
-#include <deal.II/lac/schur_complement.h>
-
-#include <deal.II/grid/grid_generator.h>
-#include <deal.II/grid/tria_boundary_lib.h>
-#include <deal.II/grid/grid_tools.h>
-
-#include <deal.II/dofs/dof_renumbering.h>
-#include <deal.II/dofs/dof_tools.h>
-
-#include <deal.II/fe/fe_q.h>
-#include <deal.II/fe/fe_system.h>
-
-#include <deal.II/numerics/vector_tools.h>
-#include <deal.II/numerics/data_out.h>
-#include <deal.II/numerics/error_estimator.h>
-
-namespace Step22
-{
-  using namespace dealii;
-
-  template <int dim>
-  class StokesProblem
-  {
-  public:
-    StokesProblem (const unsigned int degree);
-    void run ();
-
-  private:
-    void setup_dofs ();
-    void assemble_system ();
-    void solve ();
-    void get_point_value (const Point<dim> point, const int proc,
-                          Vector<double> &value) const;
-    void check_periodicity(const unsigned int cycle) const;
-    void output_results (const unsigned int refinement_cycle) const;
-    void refine_mesh ();
-
-    const unsigned int   degree;
-
-    MPI_Comm                                    mpi_communicator;
-
-    HyperShellBoundary<dim>                     boundary;
-    parallel::distributed::Triangulation<dim>   triangulation;
-    FESystem<dim>                               fe;
-    DoFHandler<dim>                             dof_handler;
-
-    ConstraintMatrix                            constraints;
-    std::vector<IndexSet>                       owned_partitioning;
-    std::vector<IndexSet>                       relevant_partitioning;
-
-    TrilinosWrappers::BlockSparsityPattern      sparsity_pattern;
-    TrilinosWrappers::BlockSparseMatrix         system_matrix;
-
-    TrilinosWrappers::MPI::BlockVector          solution;
-    TrilinosWrappers::MPI::BlockVector          system_rhs;
-
-    ConditionalOStream                          pcout;
-  };
-
-
-
-  template <int dim>
-  class BoundaryValues : public Function<dim>
-  {
-  public:
-    BoundaryValues () : Function<dim>(dim+1) {}
-
-    virtual double value (const Point<dim>   &p,
-                          const unsigned int  component = 0) const;
-
-    virtual void vector_value (const Point<dim> &p,
-                               Vector<double>   &value) const;
-  };
-
-
-  template <int dim>
-  double
-  BoundaryValues<dim>::value (const Point<dim>  &p,
-                              const unsigned int component) const
-  {
-    Assert (component < this->n_components,
-            ExcIndexRange (component, 0, this->n_components));
-
-    return (1-2*(component==0))*p[(component+1)%2]/(p[0]*p[0]+p[1]*p[1]);
-  }
-
-
-  template <int dim>
-  void
-  BoundaryValues<dim>::vector_value (const Point<dim> &p,
-                                     Vector<double>   &values) const
-  {
-    for (unsigned int c=0; c<this->n_components; ++c)
-      values(c) = BoundaryValues<dim>::value (p, c);
-  }
-
-
-
-  template <int dim>
-  class RightHandSide : public Function<dim>
-  {
-  public:
-    RightHandSide () : Function<dim>(dim+1) {}
-
-    virtual double value (const Point<dim>   &p,
-                          const unsigned int  component = 0) const;
-
-    virtual void vector_value (const Point<dim> &p,
-                               Vector<double>   &value) const;
-
-  };
-
-
-  template <int dim>
-  double
-  RightHandSide<dim>::value (const Point<dim>  &/*p*/,
-                             const unsigned int /*component*/) const
-  {
-    return 0;
-  }
-
-
-  template <int dim>
-  void
-  RightHandSide<dim>::vector_value (const Point<dim> &p,
-                                    Vector<double>   &values) const
-  {
-    for (unsigned int c=0; c<this->n_components; ++c)
-      values(c) = RightHandSide<dim>::value (p, c);
-  }
-
-
-
-
-
-  template <class Matrix, class Preconditioner>
-  class InverseMatrix : public Preconditioner
-  {
-  public:
-    InverseMatrix (const Matrix         &m,
-                   const Preconditioner &preconditioner,
-                   const IndexSet       &locally_owned,
-                   const MPI_Comm       &mpi_communicator);
-
-    void vmult (TrilinosWrappers::MPI::Vector       &dst,
-                const TrilinosWrappers::MPI::Vector &src) const;
-
-  private:
-    const SmartPointer<const Matrix> matrix;
-    const SmartPointer<const Preconditioner> preconditioner;
-
-    const MPI_Comm *mpi_communicator;
-    mutable TrilinosWrappers::MPI::Vector tmp;
-  };
-
-
-  template <class Matrix, class Preconditioner>
-  InverseMatrix<Matrix,Preconditioner>::InverseMatrix
-  (const Matrix         &m,
-   const Preconditioner &preconditioner,
-   const IndexSet       &locally_owned,
-   const MPI_Comm       &mpi_communicator)
-    :
-    matrix (&m),
-    preconditioner (&preconditioner),
-    mpi_communicator (&mpi_communicator),
-    tmp(locally_owned, mpi_communicator)
-  {}
-
-
-
-  template <class Matrix, class Preconditioner>
-  void InverseMatrix<Matrix,Preconditioner>::vmult
-  (TrilinosWrappers::MPI::Vector       &dst,
-   const TrilinosWrappers::MPI::Vector &src) const
-  {
-    SolverControl solver_control (src.size(), 1e-6*src.l2_norm(), false, false);
-    TrilinosWrappers::SolverCG    cg (solver_control,
-                                      TrilinosWrappers::SolverCG::AdditionalData());
-
-    tmp = 0.;
-    cg.solve (*matrix, tmp, src, *preconditioner);
-    dst = tmp;
-  }
-
-
-
-  template <class Preconditioner>
-  class SchurComplement : public TrilinosWrappers::SparseMatrix
-  {
-  public:
-    SchurComplement ( const TrilinosWrappers::BlockSparseMatrix &system_matrix,
-                      const InverseMatrix<TrilinosWrappers::SparseMatrix,
-                      Preconditioner> &A_inverse,
-                      const IndexSet &owned_pres,
-                      const IndexSet &relevant_pres,
-                      const MPI_Comm &mpi_communicator);
-
-    void vmult (TrilinosWrappers::MPI::Vector       &dst,
-                const TrilinosWrappers::MPI::Vector &src) const;
-
-  private:
-    const SmartPointer<const TrilinosWrappers::BlockSparseMatrix> system_matrix;
-    const SmartPointer<const InverseMatrix<TrilinosWrappers::SparseMatrix,
-          Preconditioner> > A_inverse;
-    mutable TrilinosWrappers::MPI::Vector tmp1, tmp2;
-  };
-
-
-
-  template <class Preconditioner>
-  SchurComplement<Preconditioner>::
-  SchurComplement (const TrilinosWrappers::BlockSparseMatrix &system_matrix,
-                   const InverseMatrix<TrilinosWrappers::SparseMatrix,
-                   Preconditioner> &A_inverse,
-                   const IndexSet &owned_vel,
-                   const IndexSet &relevant_vel,
-                   const MPI_Comm &mpi_communicator)
-    :
-    system_matrix (&system_matrix),
-    A_inverse (&A_inverse),
-    tmp1 (owned_vel, mpi_communicator),
-    tmp2 (tmp1)
-  {}
-
-
-  template <class Preconditioner>
-  void SchurComplement<Preconditioner>::vmult
-  (TrilinosWrappers::MPI::Vector       &dst,
-   const TrilinosWrappers::MPI::Vector &src) const
-  {
-    system_matrix->block(0,1).vmult (tmp1, src);
-    A_inverse->vmult (tmp2, tmp1);
-    system_matrix->block(1,0).vmult (dst, tmp2);
-  }
-
-
-
-
-  template <int dim>
-  StokesProblem<dim>::StokesProblem (const unsigned int degree)
-    :
-    degree (degree),
-    mpi_communicator (MPI_COMM_WORLD),
-    triangulation (mpi_communicator/*,
-                   Triangulation<dim>::maximum_smoothing*/),
-    fe (FE_Q<dim>(degree+1), dim,
-        FE_Q<dim>(degree), 1),
-    dof_handler (triangulation),
-    pcout (Utilities::MPI::this_mpi_process(mpi_communicator)
-           == 0
-           ?
-           deallog.get_file_stream()
-           :
-           std::cout,
-           (Utilities::MPI::this_mpi_process(mpi_communicator)
-            == 0))
-  {}
-
-
-
-
-  template <int dim>
-  void StokesProblem<dim>::setup_dofs ()
-  {
-    dof_handler.distribute_dofs (fe);
-
-    std::vector<unsigned int> block_component (dim+1,0);
-    block_component[dim] = 1;
-    DoFRenumbering::component_wise (dof_handler, block_component);
-
-    std::vector<types::global_dof_index> dofs_per_block (2);
-    DoFTools::count_dofs_per_block (dof_handler, dofs_per_block, block_component);
-    const unsigned int n_u = dofs_per_block[0],
-                       n_p = dofs_per_block[1];
-
-    {
-      owned_partitioning.clear();
-      IndexSet locally_owned_dofs = dof_handler.locally_owned_dofs();
-      owned_partitioning.push_back(locally_owned_dofs.get_view(0, n_u));
-      owned_partitioning.push_back(locally_owned_dofs.get_view(n_u, n_u+n_p));
-
-      relevant_partitioning.clear();
-      IndexSet locally_relevant_dofs;
-      DoFTools::extract_locally_relevant_dofs (dof_handler, locally_relevant_dofs);
-      relevant_partitioning.push_back(locally_relevant_dofs.get_view(0, n_u));
-      relevant_partitioning.push_back(locally_relevant_dofs.get_view(n_u, n_u+n_p));
-
-      constraints.clear ();
-      constraints.reinit(locally_relevant_dofs);
-
-      FEValuesExtractors::Vector velocities(0);
-      FEValuesExtractors::Scalar pressure(dim);
-
-      DoFTools::make_hanging_node_constraints (dof_handler,
-                                               constraints);
-#ifdef PERIODIC
-      VectorTools::interpolate_boundary_values (dof_handler,
-                                                3,
-                                                BoundaryValues<dim>(),
-                                                constraints,
-                                                fe.component_mask(velocities));
-#endif
-      VectorTools::interpolate_boundary_values (dof_handler,
-                                                0,
-                                                BoundaryValues<dim>(),//Functions::ZeroFunction<dim>(dim+1),
-                                                constraints,
-                                                fe.component_mask(velocities));
-      VectorTools::interpolate_boundary_values (dof_handler,
-                                                1,
-                                                BoundaryValues<dim>(),//Functions::ZeroFunction<dim>(dim+1),
-                                                constraints,
-                                                fe.component_mask(velocities));
-
-#ifdef PERIODIC
-      std::vector<GridTools::PeriodicFacePair<typename DoFHandler<dim>::cell_iterator> >
-      periodicity_vector;
-
-      FullMatrix<double> matrix(dim);
-      matrix[0][1]=1.;
-      matrix[1][0]=-1.;
-
-      std::vector<unsigned int> first_vector_components;
-      first_vector_components.push_back(0);
-
-      GridTools::collect_periodic_faces(
-        dof_handler, 2, 3, 1, periodicity_vector, Tensor<1, dim>(), matrix);
-
-      DoFTools::make_periodicity_constraints<DoFHandler<dim> >(
-        periodicity_vector, constraints, fe.component_mask(velocities),
-        first_vector_components);
-#endif
-    }
-
-    constraints.close ();
-
-    {
-      TrilinosWrappers::BlockSparsityPattern bsp
-      (owned_partitioning, owned_partitioning,
-       relevant_partitioning, mpi_communicator);
-
-      DoFTools::make_sparsity_pattern
-      (dof_handler, bsp, constraints, false,
-       Utilities::MPI::this_mpi_process(mpi_communicator));
-
-      bsp.compress();
-
-      system_matrix.reinit (bsp);
-    }
-
-    system_rhs.reinit (owned_partitioning,
-                       mpi_communicator);
-    solution.reinit (owned_partitioning, relevant_partitioning,
-                     mpi_communicator);
-  }
-
-
-
-  template <int dim>
-  void StokesProblem<dim>::assemble_system ()
-  {
-    system_matrix=0.;
-    system_rhs=0.;
-
-    QGauss<dim>   quadrature_formula(degree+2);
-
-    FEValues<dim> fe_values (fe, quadrature_formula,
-                             update_values    |
-                             update_quadrature_points  |
-                             update_JxW_values |
-                             update_gradients);
-
-    const unsigned int   dofs_per_cell   = fe.dofs_per_cell;
-
-    const unsigned int   n_q_points      = quadrature_formula.size();
-
-    FullMatrix<double>   local_matrix (dofs_per_cell, dofs_per_cell);
-    Vector<double>       local_rhs (dofs_per_cell);
-
-    std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
-
-    const RightHandSide<dim>          right_hand_side;
-    std::vector<Vector<double> >      rhs_values (n_q_points,
-                                                  Vector<double>(dim+1));
-
-    const FEValuesExtractors::Vector velocities (0);
-    const FEValuesExtractors::Scalar pressure (dim);
-
-    std::vector<SymmetricTensor<2,dim> > symgrad_phi_u (dofs_per_cell);
-    std::vector<double>                  div_phi_u   (dofs_per_cell);
-    std::vector<double>                  phi_p       (dofs_per_cell);
-
-    typename DoFHandler<dim>::active_cell_iterator
-    cell = dof_handler.begin_active(),
-    endc = dof_handler.end();
-    for (; cell!=endc; ++cell)
-      if (cell->is_locally_owned())
-        {
-          fe_values.reinit (cell);
-          local_matrix = 0;
-          local_rhs = 0;
-
-          right_hand_side.vector_value_list(fe_values.get_quadrature_points(),
-                                            rhs_values);
-
-          for (unsigned int q=0; q<n_q_points; ++q)
-            {
-              for (unsigned int k=0; k<dofs_per_cell; ++k)
-                {
-                  symgrad_phi_u[k] = fe_values[velocities].symmetric_gradient (k, q);
-                  div_phi_u[k]     = fe_values[velocities].divergence (k, q);
-                  phi_p[k]         = fe_values[pressure].value (k, q);
-                }
-
-              for (unsigned int i=0; i<dofs_per_cell; ++i)
-                {
-                  for (unsigned int j=0; j<=i; ++j)
-                    {
-                      local_matrix(i,j) += (symgrad_phi_u[i] * symgrad_phi_u[j]
-                                            - div_phi_u[i] * phi_p[j]
-                                            - phi_p[i] * div_phi_u[j]
-                                            + phi_p[i] * phi_p[j])
-                                           * fe_values.JxW(q);
-
-                    }
-
-                  const unsigned int component_i =
-                    fe.system_to_component_index(i).first;
-                  local_rhs(i) += fe_values.shape_value(i,q) *
-                                  rhs_values[q](component_i) *
-                                  fe_values.JxW(q);
-                }
-            }
-
-
-          for (unsigned int i=0; i<dofs_per_cell; ++i)
-            for (unsigned int j=i+1; j<dofs_per_cell; ++j)
-              local_matrix(i,j) = local_matrix(j,i);
-
-          cell->get_dof_indices (local_dof_indices);
-          constraints.distribute_local_to_global (local_matrix, local_rhs,
-                                                  local_dof_indices,
-                                                  system_matrix, system_rhs);
-        }
-
-    system_matrix.compress (VectorOperation::add);
-    system_rhs.compress (VectorOperation::add);
-  }
-
-
-
-
-  template <int dim>
-  void StokesProblem<dim>::solve ()
-  {
-    // Compute a reference solution
-    TrilinosWrappers::MPI::BlockVector solution_reference (owned_partitioning,
-                                                           relevant_partitioning,
-                                                           mpi_communicator);
-    {
-      deallog.push("ManualSchur");
-
-      TrilinosWrappers::PreconditionJacobi A_preconditioner;
-      A_preconditioner.initialize(system_matrix.block(0,0));
-
-      const InverseMatrix<TrilinosWrappers::SparseMatrix,
-            TrilinosWrappers::PreconditionJacobi>
-            A_inverse (system_matrix.block(0,0),
-                       A_preconditioner,
-                       owned_partitioning[0],
-                       mpi_communicator);
-
-      TrilinosWrappers::MPI::BlockVector tmp (owned_partitioning,
-                                              mpi_communicator);
-
-      {
-        TrilinosWrappers::MPI::Vector schur_rhs (owned_partitioning[1],
-                                                 mpi_communicator);
-        A_inverse.vmult (tmp.block(0), system_rhs.block(0));
-        system_matrix.block(1,0).vmult (schur_rhs, tmp.block(0));
-        schur_rhs -= system_rhs.block(1);
-
-        SchurComplement<TrilinosWrappers::PreconditionJacobi>
-        schur_complement (system_matrix, A_inverse,
-                          owned_partitioning[0], relevant_partitioning[0],
-                          mpi_communicator);
-
-        SolverControl solver_control (solution_reference.block(1).size(),
-                                      1e-6*schur_rhs.l2_norm());
-        //       TrilinosWrappers::SolverCG    cg (solver_control,
-        //                                         TrilinosWrappers::SolverCG::AdditionalData(true));
-        SolverCG<TrilinosWrappers::MPI::Vector> cg(solver_control);
-
-        TrilinosWrappers::PreconditionAMG preconditioner;
-        preconditioner.initialize (system_matrix.block(1,1));
-
-        check_solver_within_range(
-          cg.solve (schur_complement, tmp.block(1), schur_rhs, preconditioner),
-          solver_control.last_step(), 10, 12);
-
-        constraints.distribute (tmp);
-        solution_reference.block(1)=tmp.block(1);
-      }
-
-      {
-        system_matrix.block(0,1).vmult (tmp.block(0), tmp.block(1));
-        tmp.block(0) *= -1;
-        tmp.block(0) += system_rhs.block(0);
-
-        A_inverse.vmult (tmp.block(0), tmp.block(0));
-
-        constraints.distribute (tmp);
-        solution_reference.block(0)=tmp.block(0);
-      }
-
-      deallog.pop();
-    }
-
-    // Solution using SchurComplement operator
-    {
-      deallog.push("SchurComplementOp");
-
-      // Near duplicate of solver routine used in lac/schur_complement_03.cc
-      typedef TrilinosWrappers::MPI::Vector VectorType;
-      typedef TrilinosWrappers::PreconditionJacobi PreconditionerAType;
-      typedef TrilinosWrappers::PreconditionJacobi PreconditionerMType; // Solves in ~14 steps
-      // typedef TrilinosWrappers::PreconditionAMG    PreconditionerMType; // Solves in 1 step
-      typedef TrilinosWrappers::SolverCG SolverType;
-
-      // Linear operators
-      // Also note that this call to linear_operator is
-      // actually to TrilinosWrappers::linear_operator
-      const auto A  = linear_operator<VectorType>(system_matrix.block(0,0));
-      const auto B  = linear_operator<VectorType>(system_matrix.block(0,1));
-      const auto C  = linear_operator<VectorType>(system_matrix.block(1,0));
-      const auto M  = linear_operator<VectorType>(system_matrix.block(1,1)); // Mass matrix stored in this block
-      const auto D0 = null_operator<VectorType>(M);
-
-      // Cheap inverse of A
-      PreconditionerAType preconditioner_A;
-      preconditioner_A.initialize (system_matrix.block(0,0),
-                                   PreconditionerAType::AdditionalData());
-      ReductionControl solver_control_A (system_matrix.block(0,0).m(),
-                                         1e-12, 1e-6);
-      SolverType solver_A (solver_control_A);
-      const auto A_inv = inverse_operator(A, solver_A, preconditioner_A);
-
-      // Preconditioner for mass matrix
-      PreconditionerMType preconditioner_M;
-      preconditioner_M.initialize (system_matrix.block(1,1),
-                                   PreconditionerMType::AdditionalData());
-
-      // Schur complement
-      const auto S = schur_complement(A_inv,B,C,D0);
-
-      // // Proof of concept for more elaborate preconditioners constructed from
-      // // an inverse operator.
-      // // This particular implementation is rubbish, but it works. It would be
-      // // best to construct a more representative approximation S_inv_approx
-      // // and use that instead.
-      // // Approximate inverse of A
-      // typedef TrilinosWrappers::PreconditionAMG PreconditionerAExpType;
-      // PreconditionerAExpType preconditioner_A_exp;
-      // preconditioner_A_exp.initialize (system_matrix.block(0,0),
-      //                                  PreconditionerAExpType::AdditionalData());
-      // const auto A_inv_approx = linear_operator<VectorType>(system_matrix.block(0,0), preconditioner_A_exp);
-      // // Approximate Schur complement
-      // const auto S_approx = schur_complement<VectorType>(A_inv_approx,B,C,D0);
-      // IterationNumberControl solver_control_S_approx (1, 1e-12); // system_matrix.block(1,1).m()
-      // SolverType solver_S_approx (solver_control_S_approx);
-      // const auto S_inv_approx = inverse_operator(S_approx,solver_S_approx,preconditioner_M);
-
-      // Inverse of Schur complement
-      ReductionControl solver_control_S (system_matrix.block(1,1).m(),
-                                         1e-9, 1e-6);
-      SolverType solver_S (solver_control_S);
-      const auto S_inv = inverse_operator(S,solver_S,preconditioner_M);
-
-      TrilinosWrappers::MPI::BlockVector solution_tmp;
-      solution_tmp.reinit (owned_partitioning,
-                           mpi_communicator);
-      VectorType &x = solution_tmp.block(0);
-      VectorType &y = solution_tmp.block(1);
-
-      const VectorType &f = system_rhs.block(0);
-      const VectorType &g = system_rhs.block(1);
-
-      {
-        const unsigned int previous_depth = deallog.depth_file(0);
-
-        const VectorType rhs = condense_schur_rhs (A_inv,C,f,g);
-        y = S_inv * rhs;
-        x = postprocess_schur_solution (A_inv,B,y,f);
-
-        deallog.depth_file(previous_depth);
-      }
-
-
-      constraints.distribute (solution_tmp);
-      solution = solution_tmp;
-
-      deallog.pop();
-    }
-
-    // Check output of reference solution
-    // solution = solution_reference;
-
-    // The solutions will be slightly different since the solvers are
-    // configured differently. So we compute the relative error between
-    // them and expect them to remain within a reasonable margin of
-    // one another
-    // Similar goes for the output of the solution at points evaluated
-    // in check_periodicity: The solution here will not necessarily coincide
-    // exactly with the parent version of this test
-    TrilinosWrappers::MPI::BlockVector local_soln_ref(owned_partitioning,
-                                                      mpi_communicator);
-    local_soln_ref = solution_reference;
-    const double ref_local_norm = local_soln_ref.l2_norm();
-    const double ref_global_norm = std::sqrt(Utilities::MPI::sum(
-                                               ref_local_norm * ref_local_norm,
-                                               mpi_communicator));
-
-    TrilinosWrappers::MPI::BlockVector soln_diff(owned_partitioning,
-                                                 mpi_communicator);
-    soln_diff = solution;
-    soln_diff -= local_soln_ref;
-    const double local_error = soln_diff.l2_norm();
-    const double global_error = std::sqrt(Utilities::MPI::sum(
-                                            local_error * local_error,
-                                            mpi_communicator));
-
-    if (global_error/ref_global_norm < 1.e-3)
-      {
-        deallog << "Relative difference of results less than 1.e-3"
-                << std::endl;
-
-      }
-    else
-      {
-        deallog
-            << "Relative difference of results too large: "
-            << global_error/ref_global_norm
-            << std::endl;
-
-        AssertThrow(false,
-                    ExcMessage("Solution does not match reference result"));
-      }
-  }
-
-  template <int dim>
-  void StokesProblem<dim>::get_point_value
-  (const Point<dim> point, const int proc, Vector<double> &value) const
-  {
-    typename DoFHandler<dim>::active_cell_iterator cell
-      = GridTools::find_active_cell_around_point (dof_handler, point);
-
-    if (cell->is_locally_owned())
-      VectorTools::point_value (dof_handler, solution,
-                                point, value);
-
-    std::vector<double> tmp (value.size());
-    for (unsigned int i=0; i<value.size(); ++i)
-      tmp[i]=value[i];
-
-    std::vector<double> tmp2 (value.size());
-    MPI_Reduce(&(tmp[0]), &(tmp2[0]), value.size(), MPI_DOUBLE,
-               MPI_SUM, proc, mpi_communicator);
-
-    for (unsigned int i=0; i<value.size(); ++i)
-      value[i]=tmp2[i];
-  }
-
-  template <int dim>
-  void StokesProblem<dim>::check_periodicity (const unsigned int cycle) const
-  {}
-
-  template <>
-  void StokesProblem<2>::check_periodicity (const unsigned int cycle) const
-  {
-    unsigned int n_points = 4;
-    for (unsigned int i = 0; i<cycle; i++)
-      n_points*=2;
-
-    //don't test exactly at the support points, since point_value is not stable there
-    const double eps = 1./(16.*n_points);
-
-    for (unsigned int i=1; i< n_points; i++)
-      {
-        Vector<double> value1(3);
-        Vector<double> value2(3);
-
-        Point<2> point1;
-        point1(0)=0;
-        point1(1)=.5*(1.+1.*i/n_points+eps);
-        Point<2> point2;
-        point2(0)=.5*(1.+1.*i/n_points+eps);
-        point2(1)=0.;
-
-        get_point_value (point1, 0, value1);
-        get_point_value (point2, 0, value2);
-
-        if (Utilities::MPI::this_mpi_process(mpi_communicator)==0)
-          {
-            pcout << point1 << "\t" << value1[0] << "\t" << value1[1] << std::endl;
-            if (std::abs(value2[0]-value1[1])>1e-8)
-              {
-                std::cout<<point1<< "\t" << value1[1] << std::endl;
-                std::cout<<point2<< "\t" << value2[0] << std::endl;
-                Assert(false, ExcInternalError());
-              }
-            if (std::abs(value2[1]+value1[0])>1e-8)
-              {
-                std::cout<<point1<< "\t" << value1[0] << std::endl;
-                std::cout<<point2<< "\t" << value2[1] << std::endl;
-                Assert(false, ExcInternalError());
-              }
-          }
-      }
-  }
-
-
-  template <int dim>
-  void
-  StokesProblem<dim>::output_results (const unsigned int refinement_cycle)  const
-  {
-    std::vector<std::string> solution_names (dim, "velocity");
-    solution_names.push_back ("pressure");
-
-    std::vector<DataComponentInterpretation::DataComponentInterpretation>
-    data_component_interpretation
-    (dim, DataComponentInterpretation::component_is_part_of_vector);
-    data_component_interpretation
-    .push_back (DataComponentInterpretation::component_is_scalar);
-
-    DataOut<dim> data_out;
-    data_out.attach_dof_handler (dof_handler);
-    data_out.add_data_vector (solution, solution_names,
-                              DataOut<dim>::type_dof_data,
-                              data_component_interpretation);
-    Vector<float> subdomain (triangulation.n_active_cells());
-    for (unsigned int i=0; i<subdomain.size(); ++i)
-      subdomain(i) = triangulation.locally_owned_subdomain();
-    data_out.add_data_vector (subdomain, "subdomain");
-    data_out.build_patches (degree+1);
-
-    std::ostringstream filename;
-    filename << "solution-"
-             << Utilities::int_to_string (refinement_cycle, 2)
-             << "."
-             << Utilities::int_to_string (triangulation.locally_owned_subdomain(),2)
-             << ".vtu";
-
-    std::ofstream output (filename.str().c_str());
-    data_out.write_vtu (output);
-
-    if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
-      {
-        std::vector<std::string> filenames;
-        for (unsigned int i=0; i<Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD); ++i)
-          filenames.push_back (std::string("solution-") +
-                               Utilities::int_to_string (refinement_cycle, 2) +
-                               "." +
-                               Utilities::int_to_string(i, 2) +
-                               ".vtu");
-        const std::string
-        pvtu_master_filename = ("solution-" +
-                                Utilities::int_to_string (refinement_cycle, 2) +
-                                ".pvtu");
-        std::ofstream pvtu_master (pvtu_master_filename.c_str());
-        data_out.write_pvtu_record (pvtu_master, filenames);
-      }
-  }
-
-
-
-  template <int dim>
-  void
-  StokesProblem<dim>::refine_mesh ()
-  {
-
-    Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
-
-    FEValuesExtractors::Scalar pressure(dim);
-    KellyErrorEstimator<dim>::estimate (dof_handler,
-                                        QGauss<dim-1>(degree+1),
-                                        typename FunctionMap<dim>::type(),
-                                        solution,
-                                        estimated_error_per_cell,
-                                        fe.component_mask(pressure));
-
-    parallel::distributed::GridRefinement::
-    refine_and_coarsen_fixed_number (triangulation,
-                                     estimated_error_per_cell,
-                                     0.3, 0.0);
-    triangulation.execute_coarsening_and_refinement ();
-  }
-
-
-
-  template <int dim>
-  void StokesProblem<dim>::run ()
-  {
-    Point<dim> center;
-    const double inner_radius = .5;
-    const double outer_radius = 1.;
-
-    GridGenerator::quarter_hyper_shell (triangulation,
-                                        center,
-                                        inner_radius,
-                                        outer_radius,
-                                        0,
-                                        true);
-
-#ifdef PERIODIC
-    std::vector<GridTools::PeriodicFacePair<typename parallel::distributed::Triangulation<dim>::cell_iterator> >
-    periodicity_vector;
-
-    FullMatrix<double> matrix(dim);
-    matrix[0][1]=1.;
-    matrix[1][0]=-1.;
-
-    std::vector<unsigned int> first_vector_components;
-    first_vector_components.push_back(0);
-
-    GridTools::collect_periodic_faces(
-      triangulation, 2, 3, 1, periodicity_vector, Tensor<1, dim>(), matrix);
-
-    triangulation.add_periodicity(periodicity_vector);
-#endif
-
-
-    triangulation.set_boundary(0, boundary);
-    triangulation.set_boundary(1, boundary);
-
-    triangulation.refine_global (4-dim);
-
-    for (unsigned int refinement_cycle = 0; refinement_cycle<3;
-         ++refinement_cycle)
-      {
-        pcout << "Refinement cycle " << refinement_cycle << std::endl;
-
-        if (refinement_cycle > 0)
-          refine_mesh ();
-
-        setup_dofs ();
-
-        pcout << "   Assembling..." << std::endl << std::flush;
-        assemble_system ();
-
-        pcout << "   Solving..." << std::endl << std::flush;
-        solve ();
-
-#ifdef PERIODIC
-        check_periodicity(refinement_cycle);
-#endif
-        // output_results (refinement_cycle);
-
-        pcout << std::endl;
-      }
-  }
-}
-
-
-
-int main (int argc, char *argv[])
-{
-  try
-    {
-      using namespace dealii;
-      using namespace Step22;
-
-      Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1);
-
-      if (Utilities::MPI::this_mpi_process (MPI_COMM_WORLD)==0)
-        {
-          std::ofstream logfile("output");
-          deallog.attach(logfile, false);
-          {
-            StokesProblem<2> flow_problem(1);
-            flow_problem.run ();
-          }
-        }
-      else
-        {
-          StokesProblem<2> flow_problem(1);
-          flow_problem.run ();
-        }
-    }
-  catch (std::exception &exc)
-    {
-      std::cerr << std::endl << std::endl
-                << "----------------------------------------------------"
-                << std::endl;
-      std::cerr << "Exception on processing: " << std::endl
-                << exc.what() << std::endl
-                << "Aborting!" << std::endl
-                << "----------------------------------------------------"
-                << std::endl;
-
-      return 1;
-    }
-  catch (...)
-    {
-      std::cerr << std::endl << std::endl
-                << "----------------------------------------------------"
-                << std::endl;
-      std::cerr << "Unknown exception!" << std::endl
-                << "Aborting!" << std::endl
-                << "----------------------------------------------------"
-                << std::endl;
-      return 1;
-    }
-
-  return 0;
-}
diff --git a/tests/lac/schur_complement_05.with_trilinos=true.with_mpi=true.mpirun=1.output b/tests/lac/schur_complement_05.with_trilinos=true.with_mpi=true.mpirun=1.output
deleted file mode 100644 (file)
index 10fd76b..0000000
+++ /dev/null
@@ -1,43 +0,0 @@
-Refinement cycle 0
-   Assembling...
-   Solving...
-DEAL:ManualSchur::Solver stopped within 10 - 12 iterations
-DEAL::Relative difference of results less than 1.e-3
-0.00000 0.632812       -1.58049        0.00000
-0.00000 0.757812       -1.31971        0.00000
-0.00000 0.882812       -1.13281        0.00000
-
-Refinement cycle 1
-   Assembling...
-   Solving...
-DEAL:ManualSchur::Solver stopped within 10 - 12 iterations
-DEAL::Relative difference of results less than 1.e-3
-0.00000 0.566406       -1.76557        0.00000
-0.00000 0.628906       -1.59010        0.00000
-0.00000 0.691406       -1.44635        0.00000
-0.00000 0.753906       -1.32649        0.00000
-0.00000 0.816406       -1.22485        0.00000
-0.00000 0.878906       -1.13782        0.00000
-0.00000 0.941406       -1.06222        0.00000
-
-Refinement cycle 2
-   Assembling...
-   Solving...
-DEAL:ManualSchur::Solver stopped within 10 - 12 iterations
-DEAL::Relative difference of results less than 1.e-3
-0.00000 0.533203       -1.87543        0.00000
-0.00000 0.564453       -1.77163        0.00000
-0.00000 0.595703       -1.67869        0.00000
-0.00000 0.626953       -1.59502        0.00000
-0.00000 0.658203       -1.51929        0.00000
-0.00000 0.689453       -1.45043        0.00000
-0.00000 0.720703       -1.38754        0.00000
-0.00000 0.751953       -1.32988        0.00000
-0.00000 0.783203       -1.27680        0.00000
-0.00000 0.814453       -1.22782        0.00000
-0.00000 0.845703       -1.18244        0.00000
-0.00000 0.876953       -1.14033        0.00000
-0.00000 0.908203       -1.10120        0.00000
-0.00000 0.939453       -1.06444        0.00000
-0.00000 0.970703       -1.03006        0.00000
-
diff --git a/tests/lac/schur_complement_05.with_trilinos=true.with_mpi=true.mpirun=3.output b/tests/lac/schur_complement_05.with_trilinos=true.with_mpi=true.mpirun=3.output
deleted file mode 100644 (file)
index 10fd76b..0000000
+++ /dev/null
@@ -1,43 +0,0 @@
-Refinement cycle 0
-   Assembling...
-   Solving...
-DEAL:ManualSchur::Solver stopped within 10 - 12 iterations
-DEAL::Relative difference of results less than 1.e-3
-0.00000 0.632812       -1.58049        0.00000
-0.00000 0.757812       -1.31971        0.00000
-0.00000 0.882812       -1.13281        0.00000
-
-Refinement cycle 1
-   Assembling...
-   Solving...
-DEAL:ManualSchur::Solver stopped within 10 - 12 iterations
-DEAL::Relative difference of results less than 1.e-3
-0.00000 0.566406       -1.76557        0.00000
-0.00000 0.628906       -1.59010        0.00000
-0.00000 0.691406       -1.44635        0.00000
-0.00000 0.753906       -1.32649        0.00000
-0.00000 0.816406       -1.22485        0.00000
-0.00000 0.878906       -1.13782        0.00000
-0.00000 0.941406       -1.06222        0.00000
-
-Refinement cycle 2
-   Assembling...
-   Solving...
-DEAL:ManualSchur::Solver stopped within 10 - 12 iterations
-DEAL::Relative difference of results less than 1.e-3
-0.00000 0.533203       -1.87543        0.00000
-0.00000 0.564453       -1.77163        0.00000
-0.00000 0.595703       -1.67869        0.00000
-0.00000 0.626953       -1.59502        0.00000
-0.00000 0.658203       -1.51929        0.00000
-0.00000 0.689453       -1.45043        0.00000
-0.00000 0.720703       -1.38754        0.00000
-0.00000 0.751953       -1.32988        0.00000
-0.00000 0.783203       -1.27680        0.00000
-0.00000 0.814453       -1.22782        0.00000
-0.00000 0.845703       -1.18244        0.00000
-0.00000 0.876953       -1.14033        0.00000
-0.00000 0.908203       -1.10120        0.00000
-0.00000 0.939453       -1.06444        0.00000
-0.00000 0.970703       -1.03006        0.00000
-

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.