[cig-commits] commit 1993 by buerg to /var/svn/dealii/aspect
dealii.demon at gmail.com
dealii.demon at gmail.com
Fri Oct 25 09:30:54 PDT 2013
Revision 1993
Replace Uzawa by Schur-complement. Add damping of fixed point iteration for Navier-Stokes.
U trunk/aspire/include/aspect/simulator.h
U trunk/aspire/source/simulator/assembly.cc
U trunk/aspire/source/simulator/parameters.cc
U trunk/aspire/source/simulator/solver.cc
http://www.dealii.org/websvn/revision.php?repname=Aspect+Repository&path=%2F&rev=1993&peg=1993
Diff:
Modified: trunk/aspire/include/aspect/simulator.h
===================================================================
--- trunk/aspire/include/aspect/simulator.h 2013-10-25 12:58:25 UTC (rev 1992)
+++ trunk/aspire/include/aspect/simulator.h 2013-10-25 16:30:17 UTC (rev 1993)
@@ -170,6 +170,7 @@
unsigned int n_cheap_stokes_solver_steps;
double temperature_solver_tolerance;
double composition_solver_tolerance;
+ double damping;
/**
* @}
*/
@@ -961,8 +962,8 @@
LinearAlgebra::BlockVector pressure_shape_function_integrals;
- std_cxx1x::shared_ptr<LinearAlgebra::PreconditionBase> Amg_preconditioner;
- std_cxx1x::shared_ptr<LinearAlgebra::PreconditionBase> Mp_preconditioner;
+ std_cxx1x::shared_ptr<TrilinosWrappers::PreconditionBlockwiseDirect> Amg_preconditioner;
+ std_cxx1x::shared_ptr<LinearAlgebra::PreconditionILU> Mp_preconditioner;
std_cxx1x::shared_ptr<LinearAlgebra::PreconditionBase> T_preconditioner;
//TODO: use n_compositional_field separate preconditioners
std_cxx1x::shared_ptr<LinearAlgebra::PreconditionBase> C_preconditioner;
Modified: trunk/aspire/source/simulator/assembly.cc
===================================================================
--- trunk/aspire/source/simulator/assembly.cc 2013-10-25 12:58:25 UTC (rev 1992)
+++ trunk/aspire/source/simulator/assembly.cc 2013-10-25 16:30:17 UTC (rev 1993)
@@ -817,7 +817,7 @@
? -2.0/3.0 * eta + parameters.stabilization_graddiv
: parameters.stabilization_graddiv)
* (scratch.div_phi_u[i] * scratch.div_phi_u[j])
- + scratch.phi_p[i] * scratch.phi_p[j]
+ + scratch.phi_p[i] * scratch.phi_p[j] / (eta + parameters.stabilization_graddiv)
- scratch.div_phi_u[i] * scratch.phi_p[j]
- scratch.phi_p[i] * scratch.div_phi_u[j])
*
Modified: trunk/aspire/source/simulator/parameters.cc
===================================================================
--- trunk/aspire/source/simulator/parameters.cc 2013-10-25 12:58:25 UTC (rev 1992)
+++ trunk/aspire/source/simulator/parameters.cc 2013-10-25 16:30:17 UTC (rev 1993)
@@ -184,6 +184,10 @@
"the composition system gets solved. See 'linear solver "
"tolerance' for more details.");
+ prm.declare_entry ("Damping", "1",
+ Patterns::Double (0, 1),
+ "Damping factor for velocity fixed point iteration.");
+
prm.enter_subsection ("Model settings");
{
prm.declare_entry ("Include shear heating", "true",
@@ -467,6 +471,7 @@
n_cheap_stokes_solver_steps = prm.get_integer ("Number of cheap Stokes solver steps");
temperature_solver_tolerance = prm.get_double ("Temperature solver tolerance");
composition_solver_tolerance = prm.get_double ("Composition solver tolerance");
+ damping = prm.get_double ("Damping");
prm.enter_subsection ("Mesh refinement");
{
Modified: trunk/aspire/source/simulator/solver.cc
===================================================================
--- trunk/aspire/source/simulator/solver.cc 2013-10-25 12:58:25 UTC (rev 1992)
+++ trunk/aspire/source/simulator/solver.cc 2013-10-25 16:30:17 UTC (rev 1993)
@@ -34,129 +34,8 @@
namespace internal
{
using namespace dealii;
-
+
/**
- * Implement multiplication with Stokes part of system matrix
- */
- class StokesBlock
- {
- public:
- /**
- * @brief Constructor
- *
- * @param S The entire system matrix
- */
- StokesBlock (const LinearAlgebra::BlockSparseMatrix &S)
- : system_matrix(S) {}
-
- /**
- * Matrix vector product with Stokes block.
- */
- void vmult (LinearAlgebra::BlockVector &dst,
- const LinearAlgebra::BlockVector &src) const;
-
- void Tvmult (LinearAlgebra::BlockVector &dst,
- const LinearAlgebra::BlockVector &src) const;
-
- void vmult_add (LinearAlgebra::BlockVector &dst,
- const LinearAlgebra::BlockVector &src) const;
-
- void Tvmult_add (LinearAlgebra::BlockVector &dst,
- const LinearAlgebra::BlockVector &src) const;
-
- /**
- * Compute the residual with the Stokes block. In a departure from
- * the other functions, the #b variable may actually have more than
- * two blocks so that we can put it a global system_rhs vector. The
- * other vectors need to have 2 blocks only.
- */
- double residual (TrilinosWrappers::MPI::BlockVector &dst,
- const TrilinosWrappers::MPI::BlockVector &x,
- const TrilinosWrappers::MPI::BlockVector &b) const;
-
- private:
-
- /**
- * Reference to the system matrix object.
- */
- const LinearAlgebra::BlockSparseMatrix &system_matrix;
- };
-
-
-
- void StokesBlock::vmult (LinearAlgebra::BlockVector &dst,
- const LinearAlgebra::BlockVector &src) const
- {
- Assert (src.n_blocks() == 2, ExcInternalError());
- Assert (dst.n_blocks() == 2, ExcInternalError());
-
- system_matrix.block(0,0).vmult(dst.block(0), src.block(0));
- system_matrix.block(0,1).vmult_add(dst.block(0), src.block(1));
-
- system_matrix.block(1,0).vmult(dst.block(1), src.block(0));
- system_matrix.block(1,1).vmult_add(dst.block(1), src.block(1));
- }
-
- void StokesBlock::Tvmult (LinearAlgebra::BlockVector &dst,
- const LinearAlgebra::BlockVector &src) const
- {
- Assert (src.n_blocks() == 2, ExcInternalError());
- Assert (dst.n_blocks() == 2, ExcInternalError());
-
- system_matrix.block(0,0).Tvmult(dst.block(0), src.block(0));
- system_matrix.block(1,0).Tvmult_add(dst.block(0), src.block(1));
-
- system_matrix.block(0,1).Tvmult(dst.block(1), src.block(0));
- system_matrix.block(1,1).Tvmult_add(dst.block(1), src.block(1));
- }
-
- void StokesBlock::vmult_add (LinearAlgebra::BlockVector &dst,
- const LinearAlgebra::BlockVector &src) const
- {
- Assert (src.n_blocks() == 2, ExcInternalError());
- Assert (dst.n_blocks() == 2, ExcInternalError());
-
- system_matrix.block(0,0).vmult_add(dst.block(0), src.block(0));
- system_matrix.block(0,1).vmult_add(dst.block(0), src.block(1));
-
- system_matrix.block(1,0).vmult_add(dst.block(1), src.block(0));
- system_matrix.block(1,1).vmult_add(dst.block(1), src.block(1));
- }
-
- void StokesBlock::Tvmult_add (LinearAlgebra::BlockVector &dst,
- const LinearAlgebra::BlockVector &src) const
- {
- Assert (src.n_blocks() == 2, ExcInternalError());
- Assert (dst.n_blocks() == 2, ExcInternalError());
-
- system_matrix.block(0,0).Tvmult_add(dst.block(0), src.block(0));
- system_matrix.block(1,0).Tvmult_add(dst.block(0), src.block(1));
-
- system_matrix.block(0,1).Tvmult_add(dst.block(1), src.block(0));
- system_matrix.block(1,1).Tvmult_add(dst.block(1), src.block(1));
- }
-
-
-
- double StokesBlock::residual (TrilinosWrappers::MPI::BlockVector &dst,
- const TrilinosWrappers::MPI::BlockVector &x,
- const TrilinosWrappers::MPI::BlockVector &b) const
- {
- Assert (x.n_blocks() == 2, ExcInternalError());
- Assert (dst.n_blocks() == 2, ExcInternalError());
- // compute b-Ax where A is only the top left 2x2 block
- this->vmult (dst, x);
- dst.block(0).sadd (-1, 1, b.block(0));
- dst.block(1).sadd (-1, 1, b.block(1));
- // clear blocks we didn't want to fill
- for (unsigned int b=2; b<dst.n_blocks(); ++b)
- dst.block(b) = 0;
-
- return dst.l2_norm();
- }
-
-
- /**
* Implement the block Schur preconditioner for the Stokes system.
*/
template <class PreconditionerA, class PreconditionerMp>
@@ -174,8 +53,8 @@
* the matrix $A$, or only apply one preconditioner step with it.
**/
BlockSchurPreconditioner (const LinearAlgebra::BlockSparseMatrix &S,
- const PreconditionerMp *Mppreconditioner,
- const PreconditionerA *Apreconditioner,
+ const PreconditionerMp &Mppreconditioner,
+ const PreconditionerA &Apreconditioner,
const bool do_solve_A);
/**
@@ -194,8 +73,8 @@
* References to the various matrix object this preconditioner works on.
*/
const LinearAlgebra::BlockSparseMatrix &stokes_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
@@ -208,8 +87,8 @@
template <class PreconditionerA, class PreconditionerMp>
BlockSchurPreconditioner<PreconditionerA, PreconditionerMp>::
BlockSchurPreconditioner (const LinearAlgebra::BlockSparseMatrix &S,
- const PreconditionerMp *Mppreconditioner,
- const PreconditionerA *Apreconditioner,
+ const PreconditionerMp &Mppreconditioner,
+ const PreconditionerA &Apreconditioner,
const bool do_solve_A)
:
iters_A (0),
@@ -242,7 +121,7 @@
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);
+ mp_preconditioner);
dst.block(1) *= -1.0;
}
@@ -257,12 +136,13 @@
{
SolverControl solver_control(5000, utmp.l2_norm()*1e-2);
TrilinosWrappers::SolverGMRES solver(solver_control);
- solver.solve(stokes_matrix.block(0,0), dst.block(0), utmp,
- *a_preconditioner);
+ solver.solve (stokes_matrix.block (0,0), dst.block (0), utmp,
+ a_preconditioner);
iters_A += solver_control.last_step();
}
+
else
- a_preconditioner->vmult (dst.block(0), utmp);
+ a_preconditioner.vmult (dst.block(0), utmp);
}
}
@@ -370,65 +250,71 @@
distributed_stokes_solution.block (0) = current_linearization_point_velocity.block (0);
distributed_stokes_solution.block (1) = current_linearization_point_velocity.block (1);
+
+ const double solver_tolerance = std::max (parameters.linear_stokes_solver_tolerance *
+ system_rhs_velocity.l2_norm (), 1e-12);
+ PrimitiveVectorMemory<LinearAlgebra::BlockVector> mem;
+ SolverControl solver_control_cheap (parameters.n_cheap_stokes_solver_steps, solver_tolerance);
+ SolverControl solver_control_expensive (system_matrix_velocity.block (0, 1). m ()
+ + system_matrix_velocity.block (1, 0). m (), solver_tolerance);
- // extract Stokes parts of rhs vector
- LinearAlgebra::BlockVector distributed_stokes_rhs(2);
- distributed_stokes_rhs.block(0).reinit (introspection.index_sets.velocity_partitioning[0], mpi_communicator);
- distributed_stokes_rhs.block(1).reinit (introspection.index_sets.velocity_partitioning[1], mpi_communicator);
- distributed_stokes_rhs.collect_sizes();
- distributed_stokes_rhs.block (0).equ (-1.0, system_rhs_velocity.block (0));
+ try
+ {
+ // if this cheaper solver is not desired, the simply short-cut
+ // the attempt at solving with the cheaper preconditioner
+ if (parameters.n_cheap_stokes_solver_steps == 0)
+ throw SolverControl::NoConvergence(0,0);
- const double solver_tolerance = std::max (parameters.linear_stokes_solver_tolerance *
- distributed_stokes_rhs.l2_norm(), 1e-12);
- SolverControl solver_control_velocity (distributed_stokes_rhs.block (0).size (), solver_tolerance);
-
- {
- system_matrix_velocity.block (0, 1).vmult_add (distributed_stokes_rhs.block (0), distributed_stokes_solution.block (1));
- distributed_stokes_rhs.block (0) *= -1.0;
+ // otherwise give it a try with a preconditioner that consists
+ // of only a single V-cycle
+ const internal::BlockSchurPreconditioner<TrilinosWrappers::PreconditionBlockwiseDirect,
+ LinearAlgebra::PreconditionILU>
+ preconditioner (system_matrix_velocity, *Mp_preconditioner, *Amg_preconditioner, false);
+ SolverFGMRES<LinearAlgebra::BlockVector> solver (solver_control_cheap, mem,
+ SolverFGMRES<LinearAlgebra::BlockVector>::AdditionalData (30, true));
+
+ solver.solve (system_matrix_velocity, distributed_stokes_solution, system_rhs_velocity, preconditioner);
+ }
+
+ // step 1b: take the stronger solver in case
+ // the simple solver failed
+ catch (SolverControl::NoConvergence)
+ {
+ const internal::BlockSchurPreconditioner<TrilinosWrappers::PreconditionBlockwiseDirect,
+ LinearAlgebra::PreconditionILU>
+ preconditioner (system_matrix_velocity, *Mp_preconditioner, *Amg_preconditioner, true);
+ SolverFGMRES<LinearAlgebra::BlockVector> solver (solver_control_expensive, mem,
+ SolverFGMRES<LinearAlgebra::BlockVector>::AdditionalData (50, true));
- TrilinosWrappers::SolverGMRES gmres (solver_control_velocity, TrilinosWrappers::SolverGMRES::AdditionalData (50, true));
-
- gmres.solve (system_matrix_velocity.block (0, 0), distributed_stokes_solution.block (0), distributed_stokes_rhs.block (0), *Amg_preconditioner);
- }
-
- SolverControl solver_control_pressure (distributed_stokes_rhs.block (1).size (), solver_tolerance);
-
- {
- system_matrix_velocity.block (1, 0).residual (distributed_stokes_rhs.block (1), distributed_stokes_solution.block (0), system_rhs_velocity.block (1));
- distributed_stokes_rhs.block (1) *= -1e-9;
- system_matrix_velocity.block (1, 1).vmult_add (distributed_stokes_rhs.block (1), distributed_stokes_solution.block (1));
-
- TrilinosWrappers::SolverCG cg (solver_control_pressure);
-
- cg.solve (system_matrix_velocity.block (1, 1), distributed_stokes_solution.block (1), distributed_stokes_rhs.block (1), *Mp_preconditioner);
- }
-
+ solver.solve (system_matrix_velocity, distributed_stokes_solution, system_rhs_velocity, preconditioner);
+ }
// distribute hanging node and
// other constraints
current_constraints_velocity.distribute (distributed_stokes_solution);
// then copy back the solution from the temporary (non-ghosted) vector
// into the ghosted one with all solution components
- solution_velocity.block (0) = distributed_stokes_solution.block (0);
- solution_velocity.block (1) = distributed_stokes_solution.block (1);
+ solution_velocity.block (0).equ (1.0 - parameters.damping, old_solution_velocity.block (0), parameters.damping, distributed_stokes_solution.block (0));
+ solution_velocity.block (1).equ (1.0 - parameters.damping, old_solution_velocity.block (1), parameters.damping, distributed_stokes_solution.block (1));
// now rescale the pressure back to real physical units
normalize_pressure (solution_velocity);
- distributed_stokes_solution.block (0) -= current_linearization_point_velocity.block (0);
+ distributed_stokes_solution.block (0) -= old_solution_velocity.block (0);
const double nonlinear_error = distributed_stokes_solution.block (0).l2_norm ();
// print the number of iterations to screen and record it in the
// statistics file
- pcout << solver_control_velocity.last_step () << '+'
- << solver_control_pressure.last_step ()
- << std::endl;
-
- statistics.add_value("Iterations for Stokes solver",
- solver_control_velocity.last_step () + solver_control_pressure.last_step ());
-
- computing_timer.exit_section();
-
+ if (solver_control_expensive.last_step () == 0)
+ pcout << solver_control_cheap.last_step () << " iterations.";
+ else
+ pcout << solver_control_cheap.last_step () << '+'
+ << solver_control_expensive.last_step () << " iterations.";
+
+ pcout << std::endl;
+ statistics.add_value ("Iterations for Stokes solver",
+ solver_control_cheap.last_step () + solver_control_expensive.last_step ());
+ computing_timer.exit_section ();
return nonlinear_error;
}
More information about the CIG-COMMITS
mailing list