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

[dealii.demon at gmail.com]

Revision 1864

Due to lack of efficient preconditioner, use Uzawa to solve Navier-Stokes.

U   trunk/aspire/source/simulator/solver.cc


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

Diff:
Modified: trunk/aspire/source/simulator/solver.cc
===================================================================
--- trunk/aspire/source/simulator/solver.cc	2013-08-27 15:33:15 UTC (rev 1863)
+++ trunk/aspire/source/simulator/solver.cc	2013-08-30 16:47:08 UTC (rev 1864)
@@ -383,62 +383,33 @@
     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) = system_rhs_velocity.block(0);
-    distributed_stokes_rhs.block(1) = system_rhs_velocity.block(1);
+    distributed_stokes_rhs.block (0).equ (-1.0, system_rhs_velocity.block (0));
 
-    PrimitiveVectorMemory< LinearAlgebra::BlockVector > mem;
-
-    // step 1a: try if the simple and fast solver
-    // succeeds in 30 steps or less (or whatever the chosen value for the
-    // corresponding parameter is).
     const double solver_tolerance = std::max (parameters.linear_solver_tolerance *
                                               distributed_stokes_rhs.l2_norm(), 1e-12);
-    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);
-    unsigned int n_inner_iterations_A;
-
-    try
-      {
-        // if this cheaper solver is not desired
-        if (parameters.n_cheap_stokes_solver_steps == 0)
-          throw SolverControl::NoConvergence(0,0);
-        
-        // 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, &(*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,
-                      distributed_stokes_rhs, preconditioner);
-      }
-
-    // step 1b: take the stronger solver in case
-    // the simple solver failed
-    catch (SolverControl::NoConvergence)
-      {
-        const internal::BlockSchurPreconditioner<LinearAlgebra::PreconditionBase,
-                                                 LinearAlgebra::PreconditionBase>
-          preconditioner (system_matrix_velocity, &(*Mp_preconditioner), &(*Amg_preconditioner),
-                          true);
-
-        SolverFGMRES<LinearAlgebra::BlockVector>
-        solver(solver_control_expensive, mem,
-               SolverFGMRES<LinearAlgebra::BlockVector>::
-               AdditionalData(50, true));
-        solver.solve (system_matrix_velocity, distributed_stokes_solution,
-                      distributed_stokes_rhs, preconditioner);
-
-        n_inner_iterations_A = preconditioner.iters_A;
-      }
-
+    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;
+      
+      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);
+    }
+    
     // distribute hanging node and
     // other constraints
     current_constraints_velocity.distribute (distributed_stokes_solution);
@@ -449,7 +420,6 @@
     solution_velocity.block (1) = distributed_stokes_solution.block (1);
 
     // now rescale the pressure back to real physical units
-    solution_velocity.block (1) *= pressure_scaling;
     normalize_pressure (solution_velocity);
 
     distributed_stokes_solution.block (0) -= current_linearization_point_velocity.block (0);
@@ -457,17 +427,12 @@
     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
-    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 (# inner its for A: "
-            << n_inner_iterations_A  << ").";
-    pcout << std::endl;
+    pcout << solver_control_velocity.last_step () << '+'
+          << solver_control_pressure.last_step ()
+          << std::endl;
 
     statistics.add_value("Iterations for Stokes solver",
-                         solver_control_cheap.last_step() + solver_control_expensive.last_step());
+                         solver_control_velocity.last_step () + solver_control_pressure.last_step ());
 
     computing_timer.exit_section();