[cig-commits] r1385 - in branches/s-wang: for_deal.II/include/deal.II/lac for_deal.II/source/lac include/aspect source/simulator

s-wang at dealii.org s-wang at dealii.org
Mon Nov 26 07:30:30 PST 2012


Author: s-wang
Date: 2012-11-26 08:29:50 -0700 (Mon, 26 Nov 2012)
New Revision: 1385

Modified:
   branches/s-wang/for_deal.II/include/deal.II/lac/petsc_matrix_base.h
   branches/s-wang/for_deal.II/source/lac/petsc_matrix_base.cc
   branches/s-wang/include/aspect/global.h
   branches/s-wang/source/simulator/assembly.cc
   branches/s-wang/source/simulator/core.cc
   branches/s-wang/source/simulator/solver.cc
Log:
KSPSetNormType(solver_data->ksp, KSP_NORM_UNPRECONDITIONED) bas been added to the petsc solver so that it has the same norm as that of trilinos for the convergece test.

Modified: branches/s-wang/for_deal.II/include/deal.II/lac/petsc_matrix_base.h
===================================================================
--- branches/s-wang/for_deal.II/include/deal.II/lac/petsc_matrix_base.h	2012-11-21 15:10:01 UTC (rev 1384)
+++ branches/s-wang/for_deal.II/include/deal.II/lac/petsc_matrix_base.h	2012-11-26 15:29:50 UTC (rev 1385)
@@ -1142,7 +1142,7 @@
                                          * this function simply writes non-zero
                                          * elements of a matrix to the terminal.
                                          */
-      void write_ascii ();
+      void write_ascii () const;
 
                                        /**
                                         *  Returns the number bytes consumed

Modified: branches/s-wang/for_deal.II/source/lac/petsc_matrix_base.cc
===================================================================
--- branches/s-wang/for_deal.II/source/lac/petsc_matrix_base.cc	2012-11-21 15:10:01 UTC (rev 1384)
+++ branches/s-wang/for_deal.II/source/lac/petsc_matrix_base.cc	2012-11-26 15:29:50 UTC (rev 1385)
@@ -595,7 +595,7 @@
   }
 
   void
-  MatrixBase::write_ascii ()
+  MatrixBase::write_ascii () const
   {
                                        // First flush PETSc caches
     compress ();

Modified: branches/s-wang/include/aspect/global.h
===================================================================
--- branches/s-wang/include/aspect/global.h	2012-11-21 15:10:01 UTC (rev 1384)
+++ branches/s-wang/include/aspect/global.h	2012-11-26 15:29:50 UTC (rev 1385)
@@ -94,7 +94,7 @@
      */
     typedef PETScWrappers::MPI::BlockSparseMatrix BlockSparseMatrix;
 
-    typedef PETScWrappers::SolverCG SolverCG;
+//    typedef PETScWrappers::SolverCG SolverCG;
 
     /**
      * Typedef for the AMG preconditioner type used for the
@@ -112,7 +112,7 @@
      * Typedef for the Incomplete LU decomposition preconditioner used
      * for other blocks of the system matrix.
      */
-    typedef PETScWrappers::PreconditionJacobi PreconditionILU;
+    typedef PETScWrappers::PreconditionILU PreconditionILU;
   }
 }
 

Modified: branches/s-wang/source/simulator/assembly.cc
===================================================================
--- branches/s-wang/source/simulator/assembly.cc	2012-11-21 15:10:01 UTC (rev 1384)
+++ branches/s-wang/source/simulator/assembly.cc	2012-11-26 15:29:50 UTC (rev 1385)
@@ -1108,6 +1108,12 @@
     DoFTools::extract_constant_modes (dof_handler, velocity_components,
                                       constant_modes);
 
+//    for(unsigned int i=0; i<constant_modes[0].size(); i++)
+//    	std::cout << i << ": " << constant_modes[0][i] << ", " << constant_modes[1][i]<< std::endl;
+
+//    for(unsigned int i=0; i<constant_modes[1].size(); i++)
+//    	std::cout << i << ": " << constant_modes[1][i] << std::endl;
+
     Mp_preconditioner.reset (new LinearAlgebra::PreconditionILU());
     Amg_preconditioner.reset (new LinearAlgebra::PreconditionAMG());
 
@@ -1122,6 +1128,9 @@
     Amg_preconditioner->initialize (system_preconditioner_matrix.block(0,0),
                                     Amg_data);
 
+//    system_preconditioner_matrix.block(1,1).write_ascii();
+//    exit(0);
+
     rebuild_stokes_preconditioner = false;
 
     pcout << std::endl;

Modified: branches/s-wang/source/simulator/core.cc
===================================================================
--- branches/s-wang/source/simulator/core.cc	2012-11-21 15:10:01 UTC (rev 1384)
+++ branches/s-wang/source/simulator/core.cc	2012-11-26 15:29:50 UTC (rev 1385)
@@ -453,12 +453,12 @@
     // first produce some output for the screen to show where we are
     if (parameters.convert_to_years == true)
       pcout << "*** Timestep " << timestep_number
-            << ":  t=" << time/year_in_seconds
+            << ":  t=" << std::scientific << time/year_in_seconds
             << " years"
             << std::endl;
     else
       pcout << "*** Timestep " << timestep_number
-            << ":  t=" << time
+            << ":  t=" << std::scientific << time
             << " seconds"
             << std::endl;
 
@@ -1431,7 +1431,7 @@
           }
 
         postprocess ();
-        exit(0);
+//        exit(0);
 
         // see if this is a time step where additional refinement is requested
         // if so, then loop over as many times as this is necessary

Modified: branches/s-wang/source/simulator/solver.cc
===================================================================
--- branches/s-wang/source/simulator/solver.cc	2012-11-21 15:10:01 UTC (rev 1384)
+++ branches/s-wang/source/simulator/solver.cc	2012-11-26 15:29:50 UTC (rev 1385)
@@ -24,6 +24,8 @@
 #include <aspect/global.h>
 
 #include <deal.II/lac/solver_gmres.h>
+#include <deal.II/lac/solver_bicgstab.h>
+#include <deal.II/lac/solver_cg.h>
 #include <deal.II/lac/constraint_matrix.h>
 #include <deal.II/lac/trilinos_solver.h>
 #include <deal.II/lac/pointer_matrix.h>
@@ -232,8 +234,10 @@
 
       {
         SolverControl solver_control(5000, 1e-6 * src.block(1).l2_norm());
+//    	  SolverControl solver_control(5000, 1e-10);
 
-        LinearAlgebra::SolverCG solver(solver_control);
+//        PETScWrappers::SolverCG solver(solver_control);
+    	  PETScWrappers::SolverCG solver(solver_control);
 
         // Trilinos reports a breakdown
         // in case src=dst=0, even
@@ -242,9 +246,18 @@
         // iterating. We simply skip
         // solving in this case.
         if (src.block(1).l2_norm() > 1e-50 || dst.block(1).l2_norm() > 1e-50)
+        {
+//        	stokes_preconditioner_matrix.block(1,1).write_ascii();
+//        	if(dst.block(1).l2_norm()<1e-20)
+//        		dst.block(1) = 1e-10;
+//        	dst.block(1).print(std::cout,7,true,false);
+//        	src.block(1).print(std::cout,7,true,false);
           solver.solve(stokes_preconditioner_matrix.block(1,1),
                        dst.block(1), src.block(1),
                        mp_preconditioner);
+//          dst.block(1).print(std::cout,7,true,false);
+//          exit(0);
+        }
 
         dst.block(1) *= -1.0;
       }
@@ -258,7 +271,7 @@
       if (do_solve_A == true)
         {
           SolverControl solver_control(5000, utmp.l2_norm()*1e-2);
-          LinearAlgebra::SolverCG solver(solver_control);
+          PETScWrappers::SolverCG solver(solver_control);
           solver.solve(stokes_matrix.block(0,0), dst.block(0), utmp,
                        a_preconditioner);
         }
@@ -398,7 +411,7 @@
     const double solver_tolerance = std::max (parameters.linear_solver_tolerance *
                                               distributed_stokes_rhs.l2_norm(),
                                               1e-12 * initial_residual);
-    SolverControl solver_control_cheap (3000, solver_tolerance);
+    SolverControl solver_control_cheap (30, solver_tolerance);
     SolverControl solver_control_expensive (system_matrix.block(0,1).m() +
                                             system_matrix.block(1,0).m(), solver_tolerance);
 
@@ -407,9 +420,13 @@
 //    system_matrix.block(1,0).write_ascii();
 //    system_matrix.block(2,2).write_ascii();
 //    distributed_stokes_rhs.print(std::cout,7,true,false);
-//    std::cout << "initial_residual = " << initial_residual << std::endl;
-//    std::cout << "solver_tolerance = " << solver_tolerance << std::endl;
+//    system_preconditioner_matrix.block(0,0).write_ascii();
+//    system_preconditioner_matrix.block(0,1).write_ascii();
+//    system_preconditioner_matrix.block(1,0).write_ascii();
+//    system_preconditioner_matrix.block(2,2).write_ascii();
 //    exit(0);
+//    pcout << "initial_residual = " << initial_residual << std::endl;
+//    pcout << "solver_tolerance = " << solver_tolerance << std::endl;
 
     try
       {
@@ -419,10 +436,10 @@
                               *Mp_preconditioner, *Amg_preconditioner,
                               false);
 
-        SolverGMRES<LinearAlgebra::BlockVector>
+        SolverFGMRES<LinearAlgebra::BlockVector>
         solver(solver_control_cheap, mem,
-        		SolverGMRES<LinearAlgebra::BlockVector>::
-               AdditionalData(3000, true));
+        		SolverFGMRES<LinearAlgebra::BlockVector>::
+               AdditionalData(30, true));
         solver.solve(stokes_block, distributed_stokes_solution,
                      distributed_stokes_rhs, preconditioner);
       }
@@ -437,10 +454,10 @@
                               *Mp_preconditioner, *Amg_preconditioner,
                               true);
 
-        SolverGMRES<LinearAlgebra::BlockVector>
+        SolverFGMRES<LinearAlgebra::BlockVector>
         solver(solver_control_expensive, mem,
-        		SolverGMRES<LinearAlgebra::BlockVector>::
-               AdditionalData(500, true));
+        		SolverFGMRES<LinearAlgebra::BlockVector>::
+               AdditionalData(50, true));
         solver.solve(stokes_block, distributed_stokes_solution,
                      distributed_stokes_rhs, preconditioner);
       }
@@ -465,7 +482,7 @@
     // 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.";
+      pcout << solver_control_cheap.last_step()  << " iterations." << std::endl;
     else
       pcout << solver_control_cheap.last_step() << '+'
             << solver_control_expensive.last_step() << " iterations.";



More information about the CIG-COMMITS mailing list