]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Implement a few more improvements.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 25 Oct 2006 19:55:54 +0000 (19:55 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 25 Oct 2006 19:55:54 +0000 (19:55 +0000)
git-svn-id: https://svn.dealii.org/trunk@14092 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-21/Makefile
deal.II/examples/step-21/doc/intro.dox
deal.II/examples/step-21/step-21.cc

index 4aee8cd9ccb914796c707ce01b04b32f8fc612ad..26266ada331a8bbf72bb941fd6197127304594e0 100644 (file)
@@ -14,7 +14,7 @@ target = $(basename $(shell echo step-*.cc))
 # 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 = on
+debug-mode = off
 
 
 # As third field, we need to give the path to the top-level deal.II
index 41b27efc193341132770bcacfd3d7cb5c6f807df..2737f54d97fc54b2dc0f9f5c996dfd10003b552d 100644 (file)
@@ -13,7 +13,7 @@ program (besides the somewhat strange time-dependence of step-18).
 <h2>The two phase flow problem</h2>
 
 Modeling of two phase flow in porous media is important for both
-environmental remediation and the management of petroleum
+environmental remediation and the management of petroleum and groundwater
 reservoirs. Practical situations involving two phase flow include the
 dispersal of a nonaqueous phase liquid in an aquifer, or the joint
 movement of a mixture of fluids such as oil and water in a
index 4dbe61102b48ba1834806be17ef3968127c4825c..44023a8aa98319745b08fafcc0a759d8fbdcc280 100644 (file)
@@ -52,6 +52,7 @@
                                 // previous programs:
 using namespace dealii;
 
+
                                  // @sect3{The ``TwoPhaseFlowProblem'' class template}
                                  
 
@@ -68,7 +69,9 @@ class TwoPhaseFlowProblem
     void solve ();
     void compute_errors () const;
     void output_results (const unsigned int timestep_number) const;
-
+    double get_maximal_velocity () const;
+    void project_back_saturation ();
+    
     Vector<double> evaluate_solution (const Point<dim> &point) const;
     
     const unsigned int   degree;
@@ -311,13 +314,12 @@ KInverse<dim>::value_list (const std::vector<Point<dim> > &points,
 
 
 double mobility_inverse (const double S, const double vis)
-{ 
+{
   return 1.0 /(1.0/vis * S * S + (1-S) * (1-S));
 }
 
 double f_saturation(const double S, const double vis)
 {   
-
   return S*S /( S * S +vis * (1-S) * (1-S));
 }
 
@@ -851,7 +853,8 @@ void TwoPhaseFlowProblem<dim>::solve ()
     cg.solve (schur_complement, solution.block(1), schur_rhs,
               preconditioner);
   
-    std::cout << solver_control.last_step()
+    std::cout << "   "
+             << solver_control.last_step()
               << " CG Schur complement iterations to obtain convergence for pressure."
               << std::endl;
   }
@@ -1037,23 +1040,20 @@ void TwoPhaseFlowProblem<dim>::solve ()
                
       }        
     SolverControl solver_control (system_matrix.block(2,2).m(),
-                                 1e-12*system_rhs.block(2).l2_norm());
+                                 1e-8*system_rhs.block(2).l2_norm());
     SolverCG<>   cg (solver_control);
     cg.solve (system_matrix.block(2,2), solution.block(2), system_rhs.block(2),
              PreconditionIdentity());
                
        
-    std::cout << solver_control.last_step()
+    std::cout << "   "
+             << solver_control.last_step()
               << " CG iterations to obtain convergence for saturation."
               << std::endl;            
   } 
 
    
   old_solution = solution; 
-
-   
-  
 }
                                  
                                  // @sect4{TwoPhaseFlow::compute_errors}
@@ -1177,12 +1177,66 @@ void TwoPhaseFlowProblem<dim>::output_results
   data_out.build_patches (degree+1);
   
   std::ostringstream filename;
-  filename << "solution-"<< timestep_number;
+  filename << "solution-"<< timestep_number << ".vtk";
 
   std::ofstream output (filename.str().c_str());
-  data_out.write_gnuplot (output);
+  data_out.write_vtk (output);
+}
+
 
-                                  //data_out.write_vtk (output);
+
+template <int dim>
+void
+TwoPhaseFlowProblem<dim>::project_back_saturation ()
+{
+  for (unsigned int i=0; i<solution.block(dim).size(); ++i)
+    if (solution.block(dim)(i) < 0)
+      solution.block(dim)(i) = 0;
+    else
+      if (solution.block(dim)(i) > 1)
+       solution.block(dim)(i) = 1;
+
+  for (unsigned int i=0; i<solution.n_blocks(); ++i)
+    std::cout << "    sol(" << i << ")="
+             << solution.block(i).linfty_norm ()
+             << std::endl;
+}
+
+
+
+template <int dim>
+double
+TwoPhaseFlowProblem<dim>::get_maximal_velocity () const
+{
+  QGauss<dim>   quadrature_formula(degree+2); 
+  const unsigned int   n_q_points
+    = quadrature_formula.n_quadrature_points;
+
+  FEValues<dim> fe_values (fe, quadrature_formula, 
+                          update_values);
+  std::vector<Vector<double> >      old_solution_values(n_q_points, Vector<double>(dim+2));
+  double max_velocity = 0;
+  
+  typename DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+  for (; cell!=endc; ++cell)
+    {
+      fe_values.reinit (cell);
+      fe_values.get_function_values (old_solution, old_solution_values);
+
+      for (unsigned int q=0; q<n_q_points; ++q)
+       {
+         Tensor<1,dim> velocity;
+         for (unsigned int i=0; i<dim; ++i)
+           velocity[i] = old_solution_values[q](i);      
+         
+         max_velocity = std::max (max_velocity,
+                                  velocity.norm());
+       }
+    }
+
+  return max_velocity;
 }
 
 
@@ -1194,7 +1248,7 @@ void TwoPhaseFlowProblem<dim>::output_results
 template <int dim>
 void TwoPhaseFlowProblem<dim>::run () 
 {
-  std::cout<<"Solving problem in " <<dim << " space dimensions." << std::endl;
+  std::cout << "Solving problem in " <<dim << " space dimensions." << std::endl;
   
   make_grid_and_dofs();
   
@@ -1210,16 +1264,28 @@ void TwoPhaseFlowProblem<dim>::run ()
   
   unsigned int timestep_number = 1;
   
-  for ( double time = time_step; time <=1; time+=time_step,  timestep_number++)
+  for ( double time = time_step; time <= 50; time+=time_step,  timestep_number++)
     { 
-      std::cout<< "Timestep_number = "<< timestep_number<<std::endl; 
+      std::cout << "Timestep " << timestep_number
+               << " at t=" << time
+               << ", dt=" << time_step
+               << std::endl; 
       assemble_system ();
       solve ();
+      project_back_saturation ();
+      
       output_results(timestep_number);
 
       production_time.push_back (time);
       production_rate.push_back (1.0 - vfs_out/v_out);
-      std::cout<<"production_rate="<<production_rate.back()<<std::endl;       
+      std::cout << "   production_rate="<<production_rate.back()<<std::endl;
+
+      const double max_velocity = get_maximal_velocity();
+      std::cout << "   max velocity = " << max_velocity
+               << std::endl;
+      
+//       time_step = std::pow(0.5, double(n_refinement_steps)) /
+//               max_velocity / 4;
     }
 
   std::ofstream production_history ("production_history");

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.