[cig-commits] commit: Initial commit of files from Wolfgang with deal II directory modified
Mercurial
hg at geodynamics.org
Mon Feb 21 00:15:44 PST 2011
changeset: 0:62413b0e7b38
user: Walter Landry <wlandry at caltech.edu>
date: Fri Feb 18 16:50:04 2011 -0800
files: Makefile madds-1.cc madds-1.prm
description:
Initial commit of files from Wolfgang with deal II directory modified
diff -r 000000000000 -r 62413b0e7b38 Makefile
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/Makefile Fri Feb 18 16:50:04 2011 -0800
@@ -0,0 +1,143 @@
+# $Id: Makefile 22601 2010-11-04 03:26:19Z bangerth $
+
+
+# For the small projects Makefile, you basically need to fill in only
+# four fields.
+#
+# The first is the name of the application. It is assumed that the
+# application name is the same as the base file name of the single C++
+# file from which the application is generated.
+target = $(basename $(shell echo madds-*.cc))
+
+# The second field determines whether you want to run your program in
+# debug or optimized mode. The latter is significantly faster, but no
+# run-time checking of parameters and internal states is performed, so
+# you should set this value to `on' while you develop your program,
+# and to `off' when running production computations.
+debug-mode = off
+
+
+# As third field, we need to give the path to the top-level deal.II
+# directory. You need to adjust this to your needs. Since this path is
+# probably the most often needed one in the Makefile internals, it is
+# designated by a single-character variable, since that can be
+# reference using $D only, i.e. without the parentheses that are
+# required for most other parameters, as e.g. in $(target).
+D = ../deal.II
+
+
+# The last field specifies the names of data and other files that
+# shall be deleted when calling `make clean'. Object and backup files,
+# executables and the like are removed anyway. Here, we give a list of
+# files in the various output formats that deal.II supports.
+clean-up-files = *gmv *gnuplot *gpl *eps *pov *vtk *ucd *.d2
+
+
+
+
+#
+#
+# Usually, you will not need to change anything beyond this point.
+#
+#
+# The next statement tells the `make' program where to find the
+# deal.II top level directory and to include the file with the global
+# settings
+include $D/common/Make.global_options
+
+
+# Since the whole project consists of only one file, we need not
+# consider difficult dependencies. We only have to declare the
+# libraries which we want to link to the object file. deal.II has two
+# libraries: one for the debug mode version of the
+# application and one for optimized mode.
+libs.g := $(lib-deal2.g)
+libs.o := $(lib-deal2.o)
+
+
+# We now use the variable defined above to switch between debug and
+# optimized mode to select the set of libraries to link with. Included
+# in the list of libraries is the name of the object file which we
+# will produce from the single C++ file. Note that by default we use
+# the extension .g.o for object files compiled in debug mode and .o for
+# object files in optimized mode (or whatever local default on your
+# system is instead of .o)
+ifeq ($(debug-mode),on)
+ libraries = $(target).g.$(OBJEXT) $(libs.g)
+else
+ libraries = $(target).$(OBJEXT) $(libs.o)
+endif
+
+
+# Now comes the first production rule: how to link the single object
+# file produced from the single C++ file into the executable. Since
+# this is the first rule in the Makefile, it is the one `make' selects
+# if you call it without arguments.
+$(target)$(EXEEXT) : $(libraries)
+ @echo ============================ Linking $@
+ @$(CXX) -o $@ $^ $(LIBS) $(LDFLAGS)
+
+
+# To make running the application somewhat independent of the actual
+# program name, we usually declare a rule `run' which simply runs the
+# program. You can then run it by typing `make run'. This is also
+# useful if you want to call the executable with arguments which do
+# not change frequently. You may then want to add them to the
+# following rule:
+run: $(target)$(EXEEXT)
+ @echo ============================ Running $<
+ @./$(target)$(EXEEXT)
+
+
+# As a last rule to the `make' program, we define what to do when
+# cleaning up a directory. This usually involves deleting object files
+# and other automatically created files such as the executable itself,
+# backup files, and data files. Since the latter are not usually quite
+# diverse, you needed to declare them at the top of this file.
+clean:
+ -rm -f *.$(OBJEXT) *~ Makefile.dep $(target)$(EXEEXT) $(clean-up-files)
+
+
+# Since we have not yet stated how to make an object file from a C++
+# file, we should do so now. Since the many flags passed to the
+# compiler are usually not of much interest, we suppress the actual
+# command line using the `at' sign in the first column of the rules
+# and write the string indicating what we do instead.
+./%.g.$(OBJEXT) :
+ @echo "==============debug========= $(<F) -> $@"
+ @$(CXX) $(CXXFLAGS.g) -c $< -o $@
+./%.$(OBJEXT) :
+ @echo "==============optimized===== $(<F) -> $@"
+ @$(CXX) $(CXXFLAGS.o) -c $< -o $@
+
+
+# The following statement tells make that the rules `run' and `clean'
+# are not expected to produce files of the same name as Makefile rules
+# usually do.
+.PHONY: run clean
+
+
+# Finally there is a rule which you normally need not care much about:
+# since the executable depends on some include files from the library,
+# besides the C++ application file of course, it is necessary to
+# re-generate the executable when one of the files it depends on has
+# changed. The following rule creates a dependency file
+# `Makefile.dep', which `make' uses to determine when to regenerate
+# the executable. This file is automagically remade whenever needed,
+# i.e. whenever one of the cc-/h-files changed. Make detects whether
+# to remake this file upon inclusion at the bottom of this file.
+#
+# If the creation of Makefile.dep fails, blow it away and fail
+Makefile.dep: $(target).cc Makefile \
+ $(shell echo $D/include/deal.II/*/*.h)
+ @echo ============================ Remaking $@
+ @$D/common/scripts/make_dependencies $(INCLUDE) -B. $(target).cc \
+ > $@ \
+ || (rm -f $@ ; false)
+ @if test -s $@ ; then : else rm $@ ; fi
+
+# To make the dependencies known to `make', we finally have to include
+# them:
+include Makefile.dep
+
+
diff -r 000000000000 -r 62413b0e7b38 madds-1.cc
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/madds-1.cc Fri Feb 18 16:50:04 2011 -0800
@@ -0,0 +1,828 @@
+/* $Id: madds-1.cc 2224 2011-01-13 20:24:28Z bangerth $ */
+/* Author: Wolfgang Bangerth, Texas A&M University, 2008--2011 */
+
+/* $Id: madds-1.cc 2224 2011-01-13 20:24:28Z bangerth $ */
+/* */
+/* Copyright (C) 2008, 2009, 2010, 2011 by the deal.II authors */
+/* */
+/* This file is subject to QPL and may not be distributed */
+/* without copyright and license information. Please refer */
+/* to the file deal.II/doc/license.html for the text and */
+/* further information on this license. */
+
+
+ // @sect3{Include files}
+
+ // This program being a variant of step-22,
+ // the include files are not significantly
+ // different. In general, you should be
+ // familiar with step-22 before reading
+ // through the current program.
+#include <base/quadrature_lib.h>
+#include <base/logstream.h>
+#include <base/function.h>
+#include <base/utilities.h>
+#include <base/parameter_handler.h>
+
+#include <lac/block_vector.h>
+#include <lac/full_matrix.h>
+#include <lac/block_sparse_matrix.h>
+#include <lac/solver_cg.h>
+#include <lac/solver_gmres.h>
+#include <lac/precondition.h>
+#include <lac/constraint_matrix.h>
+
+#include <grid/tria.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_boundary_lib.h>
+#include <grid/grid_tools.h>
+#include <grid/grid_refinement.h>
+
+#include <dofs/dof_handler.h>
+#include <dofs/dof_renumbering.h>
+#include <dofs/dof_accessor.h>
+#include <dofs/dof_tools.h>
+
+#include <fe/fe_q.h>
+#include <fe/fe_system.h>
+#include <fe/fe_values.h>
+#include <fe/mapping_q1.h>
+
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+#include <numerics/data_out.h>
+#include <numerics/error_estimator.h>
+
+#include <lac/sparse_direct.h>
+
+#include <lac/sparse_ilu.h>
+
+#include <fstream>
+#include <sstream>
+
+using namespace dealii;
+
+
+ // @sect3{Run-time parameters}
+
+ // This program uses the ParameterHandler
+ // class to deal with a number of run-time
+ // parameters. The mechanism is essentially
+ // the same as that already used in a number
+ // of previous programs such as step-28,
+ // step-29, step-33 or step-34. The
+ // ParameterHandler class also has a lengthy
+ // introduction into its use. That said, the
+ // following should also be mostly
+ // self-explanatory: we declare a class that
+ // contains the various run-time parameters
+ // as member variables. We then need one
+ // members function that explains to the
+ // ParameterHandler class what variables
+ // exist (and what kind of values they can
+ // take, such as real numbers or integers
+ // zero or greater), and one member function
+ // that recovers the values of these
+ // parameters from the ParameterHandler
+ // object once they have been read from an
+ // input file.
+struct Parameters
+{
+ Parameters (const std::string &filename);
+
+ double eta_jump, eta;
+ unsigned int fe_degree;
+
+ unsigned int initial_refinement, refinement_cycles;
+
+ private:
+ static void declare_parameters (ParameterHandler &prm);
+ void read_parameters (ParameterHandler &prm);
+};
+
+
+Parameters::Parameters (const std::string &filename)
+{
+ ParameterHandler prm;
+ declare_parameters (prm);
+
+ prm.read_input (filename);
+ read_parameters (prm);
+}
+
+
+void
+Parameters::declare_parameters (ParameterHandler &prm)
+{
+ prm.declare_entry ("Jump in eta", "1",
+ Patterns::Double (0),
+ "Jump in the viscosity in the interior of the domain.");
+ prm.declare_entry ("eta", "1",
+ Patterns::Double (0),
+ "Viscosity parameter in the Stokes equation.");
+ prm.declare_entry ("Finite element degree", "1",
+ Patterns::Integer (1),
+ "Polynomial degree of the pressure element; the "
+ "velocity element has a polynomial degree one higher.");
+ prm.declare_entry ("Initial refinement", "4",
+ Patterns::Integer (1),
+ "Number of initial global refinement steps.");
+ prm.declare_entry ("Refinement cycles", "4",
+ Patterns::Integer (1),
+ "Number of adaptive refinement cycles.");
+}
+
+
+void
+Parameters::read_parameters (ParameterHandler &prm)
+{
+ eta_jump = prm.get_double ("Jump in eta");
+ eta = prm.get_double ("eta");
+ fe_degree = prm.get_integer ("Finite element degree");
+ initial_refinement = prm.get_integer ("Initial refinement");
+ refinement_cycles = prm.get_integer ("Refinement cycles");
+}
+
+
+
+
+ // @sect3{Support classes for linear solvers}
+
+ // The following classes are taken verbatim
+ // from step-22, or from the discussion of
+ // further improvements in the results
+ // section of step-22. Take a look there for
+ // a more thorough discussion of what is
+ // going on here:
+template <int dim>
+struct InnerPreconditioner;
+
+template <>
+struct InnerPreconditioner<2>
+{
+ typedef SparseDirectUMFPACK type;
+};
+
+template <>
+struct InnerPreconditioner<3>
+{
+ typedef SparseILU<double> type;
+};
+
+
+
+template <class Matrix, class Preconditioner>
+class InverseMatrix : public Subscriptor
+{
+ public:
+ InverseMatrix (const Matrix &m,
+ const Preconditioner &preconditioner);
+
+ void vmult (Vector<double> &dst,
+ const Vector<double> &src) const;
+
+ private:
+ const SmartPointer<const Matrix> matrix;
+ const SmartPointer<const Preconditioner> preconditioner;
+};
+
+
+template <class Matrix, class Preconditioner>
+InverseMatrix<Matrix,Preconditioner>::InverseMatrix (const Matrix &m,
+ const Preconditioner &preconditioner)
+ :
+ matrix (&m),
+ preconditioner (&preconditioner)
+{}
+
+
+
+template <class Matrix, class Preconditioner>
+void
+InverseMatrix<Matrix,Preconditioner>::vmult (Vector<double> &dst,
+ const Vector<double> &src) const
+{
+ SolverControl solver_control (src.size(), 1e-6*src.l2_norm());
+ SolverCG<> cg (solver_control);
+
+ dst = 0;
+
+ cg.solve (*matrix, dst, src, *preconditioner);
+}
+
+
+
+template <class PreconditionerA, class PreconditionerMp>
+class BlockSchurPreconditioner : public Subscriptor
+{
+ public:
+ BlockSchurPreconditioner (const BlockSparseMatrix<double> &S,
+ const InverseMatrix<SparseMatrix<double>,PreconditionerMp> &Mpinv,
+ const PreconditionerA &Apreconditioner);
+
+ void vmult (BlockVector<double> &dst,
+ const BlockVector<double> &src) const;
+
+ private:
+ const SmartPointer<const BlockSparseMatrix<double> > system_matrix;
+ const SmartPointer<const InverseMatrix<SparseMatrix<double>,
+ PreconditionerMp > > m_inverse;
+ const PreconditionerA &a_preconditioner;
+
+ mutable Vector<double> tmp;
+};
+
+
+template <class PreconditionerA, class PreconditionerMp>
+BlockSchurPreconditioner<PreconditionerA, PreconditionerMp>::
+BlockSchurPreconditioner
+(const BlockSparseMatrix<double> &S,
+ const InverseMatrix<SparseMatrix<double>,PreconditionerMp> &Mpinv,
+ const PreconditionerA &Apreconditioner)
+ :
+ system_matrix (&S),
+ m_inverse (&Mpinv),
+ a_preconditioner (Apreconditioner),
+ tmp (S.block(1,1).m())
+{}
+
+
+template <class PreconditionerA, class PreconditionerMp>
+void
+BlockSchurPreconditioner<PreconditionerA, PreconditionerMp>::
+vmult (BlockVector<double> &dst,
+ const BlockVector<double> &src) const
+{
+ a_preconditioner.vmult (dst.block(0), src.block(0));
+ system_matrix->block(1,0).residual(tmp, dst.block(0), src.block(1));
+ tmp *= -1;
+ m_inverse->vmult (dst.block(1), tmp);
+}
+
+
+ // @sect3{Equation data}
+
+ // Next a few classes that describe the
+ // boundary values at the top and on the
+ // sides , as well as the right hand side for
+ // this problem. In the computation of
+ // boundary values at the side, keep in mind
+ // that the problem has <code>dim+1</code>
+ // components. We therefore first have to
+ // isolate the velocity components
+ // (i.e. where <code>component@<dim</code>)
+ // and then set the individual components as
+ // appropriate. Note that the function should
+ // never be evaluated for pressure values
+ // since we don't need Dirichlet values for
+ // those. The definition of the depth
+ // variable $z$ is discussed in the
+ // introduction.
+template <int dim>
+class TopBoundaryValues : public Function<dim>
+{
+ public:
+ TopBoundaryValues ()
+ :
+ 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
+TopBoundaryValues<dim>::value (const Point<dim> &p,
+ const unsigned int component) const
+{
+ Assert (component < this->n_components,
+ ExcIndexRange (component, 0, this->n_components));
+
+ if (component == 0)
+ return 1;
+ return 0;
+}
+
+
+
+template <int dim>
+void
+TopBoundaryValues<dim>::vector_value (const Point<dim> &p,
+ Vector<double> &values) const
+{
+ for (unsigned int c=0; c<this->n_components; ++c)
+ values(c) = TopBoundaryValues<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);
+}
+
+
+
+
+ // @sect3{The StokesProblem class}
+
+ // Now, finally, for the main class of this
+ // program. Again, it is pretty much verbatim
+ // lifted from step-22 and the extensions
+ // section of that program. The main
+ // structural difference is that the class
+ // now has a member variable
+ // <code>parameters</code> that we initialize
+ // in the constructor by passing to it a
+ // filename from which run-time parameters
+ // should be read. The remaining (minor)
+ // differences will be discussed with the
+ // individual member functions.
+template <int dim>
+class StokesProblem
+{
+ public:
+ StokesProblem (const std::string &prm_filename);
+ void run ();
+
+ private:
+ double get_viscosity (const Point<dim> &p) const;
+
+ void setup_dofs ();
+ void assemble_system ();
+ void solve ();
+ void output_results (const unsigned int refinement_cycle) const;
+ void refine_mesh ();
+
+ const Parameters parameters;
+
+ Triangulation<dim> triangulation;
+ FESystem<dim> fe;
+ DoFHandler<dim> dof_handler;
+
+ ConstraintMatrix constraints;
+
+ BlockSparsityPattern sparsity_pattern;
+ BlockSparseMatrix<double> system_matrix;
+
+ BlockVector<double> solution;
+ BlockVector<double> system_rhs;
+
+ std_cxx1x::shared_ptr<typename InnerPreconditioner<dim>::type> A_preconditioner;
+};
+
+
+
+template <int dim>
+StokesProblem<dim>::StokesProblem (const std::string &prm_filename)
+ :
+ parameters (prm_filename),
+ triangulation (Triangulation<dim>::maximum_smoothing),
+ fe (FE_Q<dim>(parameters.fe_degree+1), dim,
+ FE_Q<dim>(parameters.fe_degree), 1),
+ dof_handler (triangulation)
+{}
+
+
+
+ // The function that sets up degrees of
+ // freedom is the same as before, except that
+ // it now has to interpolate boundary values
+ // also on the side boundary:
+template <int dim>
+void StokesProblem<dim>::setup_dofs ()
+{
+ A_preconditioner.reset ();
+ system_matrix.clear ();
+
+ dof_handler.distribute_dofs (fe);
+ DoFRenumbering::Cuthill_McKee (dof_handler);
+
+ std::vector<unsigned int> block_component (dim+1,0);
+ block_component[dim] = 1;
+ DoFRenumbering::component_wise (dof_handler, block_component);
+
+ {
+ constraints.clear ();
+ std::vector<bool> component_mask (dim+1, true);
+ component_mask[dim] = false;
+ VectorTools::interpolate_boundary_values (dof_handler,
+ 1,
+ TopBoundaryValues<dim>(),
+ constraints,
+ component_mask);
+ VectorTools::interpolate_boundary_values (dof_handler,
+ 0,
+ ZeroFunction<dim>(dim+1),
+ constraints,
+ component_mask);
+ DoFTools::make_hanging_node_constraints (dof_handler,
+ constraints);
+ }
+
+ constraints.close ();
+
+ std::vector<unsigned int> 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];
+
+ std::cout << " Number of active cells: "
+ << triangulation.n_active_cells()
+ << std::endl
+ << " Number of degrees of freedom: "
+ << dof_handler.n_dofs()
+ << " (" << n_u << '+' << n_p << ')'
+ << std::endl;
+
+ {
+ BlockCompressedSimpleSparsityPattern csp (2,2);
+
+ csp.block(0,0).reinit (n_u, n_u);
+ csp.block(1,0).reinit (n_p, n_u);
+ csp.block(0,1).reinit (n_u, n_p);
+ csp.block(1,1).reinit (n_p, n_p);
+
+ csp.collect_sizes();
+
+ DoFTools::make_sparsity_pattern (dof_handler, csp, constraints, false);
+ sparsity_pattern.copy_from (csp);
+ }
+
+ system_matrix.reinit (sparsity_pattern);
+
+ solution.reinit (2);
+ solution.block(0).reinit (n_u);
+ solution.block(1).reinit (n_p);
+ solution.collect_sizes ();
+
+ system_rhs.reinit (2);
+ system_rhs.block(0).reinit (n_u);
+ system_rhs.block(1).reinit (n_p);
+ system_rhs.collect_sizes ();
+}
+
+
+
+template <int dim>
+double StokesProblem<dim>::get_viscosity (const Point<dim> &p) const
+{
+ bool inside_high_viscosity_region = true;
+ for (unsigned int d=0; d<dim; ++d)
+ if ((p[d] < 1./3)
+ ||
+ (p[d] > 2./3))
+ {
+ inside_high_viscosity_region = false;
+ break;
+ }
+ if (inside_high_viscosity_region)
+ return parameters.eta_jump * parameters.eta;
+ else
+ return parameters.eta;
+}
+
+
+template <int dim>
+void StokesProblem<dim>::assemble_system ()
+{
+ system_matrix=0;
+ system_rhs=0;
+
+ QGauss<dim> quadrature_formula(parameters.fe_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<unsigned int> 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> > phi_grads_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)
+ {
+ 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)
+ {
+ phi_grads_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);
+ }
+
+ const double viscosity = get_viscosity (fe_values.get_quadrature_points()[q]);
+
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ for (unsigned int j=0; j<=i; ++j)
+ local_matrix(i,j) += (2 * viscosity *
+ (phi_grads_u[i] * phi_grads_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);
+ }
+
+ std::cout << " Computing preconditioner..." << std::endl << std::flush;
+
+ A_preconditioner
+ = std_cxx1x::shared_ptr<typename InnerPreconditioner<dim>::type>(new typename InnerPreconditioner<dim>::type());
+ A_preconditioner->initialize (system_matrix.block(0,0),
+ typename InnerPreconditioner<dim>::type::AdditionalData());
+
+}
+
+
+
+ // The following is the second difference:
+ // the solver now doesn't use the Schur
+ // complement any more, as we did in
+ // step-22. Rather, we use the solver and
+ // preconditioner discussed in the results
+ // section of step-22:
+template <int dim>
+void StokesProblem<dim>::solve ()
+{
+ SparseMatrix<double> pressure_mass_matrix;
+ pressure_mass_matrix.reinit(sparsity_pattern.block(1,1));
+ pressure_mass_matrix.copy_from(system_matrix.block(1,1));
+ system_matrix.block(1,1) = 0;
+
+ SparseILU<double> pmass_preconditioner;
+ pmass_preconditioner.initialize (pressure_mass_matrix,
+ SparseILU<double>::AdditionalData());
+
+ InverseMatrix<SparseMatrix<double>,SparseILU<double> >
+ m_inverse (pressure_mass_matrix, pmass_preconditioner);
+
+ BlockSchurPreconditioner<typename InnerPreconditioner<dim>::type,
+ SparseILU<double> >
+ preconditioner (system_matrix, m_inverse, *A_preconditioner);
+
+ SolverControl solver_control (system_matrix.m(),
+ 1e-6*system_rhs.l2_norm());
+ GrowingVectorMemory<BlockVector<double> > vector_memory;
+ SolverGMRES<BlockVector<double> >::AdditionalData gmres_data;
+ gmres_data.max_n_tmp_vectors = 100;
+
+ SolverGMRES<BlockVector<double> > gmres(solver_control, vector_memory,
+ gmres_data);
+
+ gmres.solve(system_matrix, solution, system_rhs,
+ preconditioner);
+
+ constraints.distribute (solution);
+
+ std::cout << " "
+ << solver_control.last_step()
+ << " block GMRES iterations"
+ << std::endl;
+}
+
+
+
+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> viscosities (triangulation.n_active_cells());
+ {
+ unsigned int index = 0;
+ for (typename Triangulation<dim>::active_cell_iterator
+ cell = triangulation.begin_active();
+ cell != triangulation.end(); ++cell, ++index)
+ viscosities(index) = get_viscosity (cell->center());
+ }
+ data_out.add_data_vector (viscosities, "viscosity_at_cell_center");
+ data_out.build_patches ();
+
+ std::ostringstream filename;
+ filename << "solution-"
+ << Utilities::int_to_string (refinement_cycle, 2)
+ << ".vtk";
+
+ std::ofstream output (filename.str().c_str());
+ data_out.write_vtk (output);
+}
+
+
+
+template <int dim>
+void
+StokesProblem<dim>::refine_mesh ()
+{
+ Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
+
+ std::vector<bool> component_mask (dim+1, false);
+ component_mask[dim] = true;
+ KellyErrorEstimator<dim>::estimate (dof_handler,
+ QGauss<dim-1>(parameters.fe_degree+1),
+ typename FunctionMap<dim>::type(),
+ solution,
+ estimated_error_per_cell,
+ component_mask);
+
+ GridRefinement::refine_and_coarsen_fixed_number (triangulation,
+ estimated_error_per_cell,
+ 0.3, 0.0);
+ triangulation.execute_coarsening_and_refinement ();
+}
+
+
+
+ // The run function is again almost
+ // identical to the one in
+ // step-22. The mesh is set up with a
+ // slightly different geometry, but
+ // the rest is the same:
+template <int dim>
+void StokesProblem<dim>::run ()
+{
+ GridGenerator::hyper_cube (triangulation);
+
+ for (typename Triangulation<dim>::active_cell_iterator
+ cell = triangulation.begin_active();
+ cell != triangulation.end(); ++cell)
+ for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+ if (triangulation.begin_active()->face(f)->center()[dim-1] == 1)
+ cell->face(f)->set_all_boundary_indicators(1);
+ triangulation.refine_global (parameters.initial_refinement);
+
+ for (unsigned int refinement_cycle = 0;
+ refinement_cycle<parameters.refinement_cycles;
+ ++refinement_cycle)
+ {
+ std::cout << "Refinement cycle " << refinement_cycle << std::endl;
+
+ if (refinement_cycle > 0)
+ refine_mesh ();
+
+ setup_dofs ();
+
+ std::cout << " Assembling..." << std::endl << std::flush;
+ assemble_system ();
+
+ std::cout << " Solving..." << std::flush;
+ solve ();
+
+ output_results (refinement_cycle);
+
+ std::cout << std::endl;
+ }
+}
+
+
+ // @sect3{The main function}
+ //
+ // The main function follows the usual
+ // scheme. As in several other programs, if
+ // no argument is given on the command line,
+ // we use <code>madds-1.prm</code> as the
+ // input file from which to read
+ // parameters. If one argument is given we
+ // use that file for parameters, and if more
+ // than one argument is given we produce an
+ // error.
+int main (int argc, char *argv[])
+{
+ try
+ {
+ deallog.depth_console (0);
+
+ if (argc > 2)
+ {
+ std::cout << "Usage:" << argv[0] << " input_file" << std::endl;
+ std::exit(1);
+ }
+
+ StokesProblem<2> flow_problem(argc==2 ?
+ argv[1] :
+ "madds-1.prm");
+ 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 -r 000000000000 -r 62413b0e7b38 madds-1.prm
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/madds-1.prm Fri Feb 18 16:50:04 2011 -0800
@@ -0,0 +1,17 @@
+# Listing of Parameters
+# ---------------------
+# Polynomial degree of the pressure element; the velocity element has a
+# polynomial degree one higher.
+set Finite element degree = 1
+
+# Number of initial global refinement steps.
+set Initial refinement = 2
+
+# Number of adaptive refinement cycles.
+set Refinement cycles = 6
+
+# Base viscosity
+set eta = 1
+
+# Jump in the viscosity in the interior of the domain.
+set Jump in eta = 1e0
More information about the CIG-COMMITS
mailing list