From: bangerth Date: Thu, 26 Oct 2006 21:22:00 +0000 (+0000) Subject: Check the thing with the random coefficient back in. The problem was that we re-initi... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=4c0a547141811dc493929a363fe0bee7d8011c8d;p=dealii-svn.git Check the thing with the random coefficient back in. The problem was that we re-initialized the random coefficient in each iteration again and again. git-svn-id: https://svn.dealii.org/trunk@14098 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-21/step-21.cc b/deal.II/examples/step-21/step-21.cc index 81e0366cfd..439a69ef8f 100644 --- a/deal.II/examples/step-21/step-21.cc +++ b/deal.II/examples/step-21/step-21.cc @@ -1,4 +1,4 @@ - // @sect3{Include files} + // @sect3{Include files} // This program is an daptation of step-20 // and includes some technique of DG method from step-12 @@ -24,7 +24,7 @@ #include #include #include - //The Discontinuous Galerkin finite element is declared: + //The Discontinuous Galerkin finite element is declared: #include #include @@ -36,7 +36,7 @@ #include #include #include - // The Raviart-Thomas finite element is declared: + // The Raviart-Thomas finite element is declared: #include // In this program, we use a tensorial @@ -48,6 +48,8 @@ // offers such functionality: #include + // The last step is as in all + // previous programs: using namespace dealii; @@ -67,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 evaluate_solution (const Point &point) const; const unsigned int degree; @@ -83,7 +87,6 @@ class TwoPhaseFlowProblem const unsigned int n_refinement_steps; double time_step; - double epsilon; double vis; double vfs_out; double v_out; @@ -98,7 +101,7 @@ class TwoPhaseFlowProblem //{Right hand side, boundary values and initial values} - // we define the template for pressure right-hand side(source function) + // we define the template for pressure right-hand side(source function) //and boundary values for pressure and saturation // initial values for saturation. @@ -198,11 +201,11 @@ SaturationBoundaryValues::vector_value (const Point &p, template double InitialValues::value (const Point &p, - const unsigned int component) const + const unsigned int component) const { - if(component::vector_value (const Point &p, // ``KInverse''.As in introduction, ' // assume the heterogeneous is isotropic, // so it is a scalar multipy the identity matrix. - //DealII has a base class not only for + //DealII has a base class not only for // scalar and generally vector-valued // functions (the ``Function'' base // class) but also for functions that @@ -258,34 +261,35 @@ InitialValues::vector_value (const Point &p, // of the function in the second // argument, a list of tensors: template -class KInverse //: public TensorFunction<2,dim> +class KInverse : public TensorFunction<2,dim> { public: - KInverse () - {} + KInverse (); - /*virtual*/ void value_list (const std::vector > &points, - std::vector > &values, - const double epsilon) const; + virtual void value_list (const std::vector > &points, + std::vector > &values) const; + + private: + std::vector > centers; }; template -class Coefficient : public Function +KInverse::KInverse () { - public: - Coefficient () : Function(1) {}; - - virtual double value (const Point &p, - const unsigned int component = 0 ) const; -}; + const unsigned int N = 40; + centers.resize (N); + for (unsigned int i=0; i void KInverse::value_list (const std::vector > &points, - std::vector > &values, - const double epsilon) const + std::vector > &values) const { Assert (points.size() == values.size(), ExcDimensionMismatch (points.size(), values.size())); @@ -294,37 +298,31 @@ KInverse::value_list (const std::vector > &points, { values[p].clear (); - const double distance_to_flowline - = std::fabs(points[p][1]-0.5-0.1*std::sin(10*points[p][0])); + double permeability = 0; + for (unsigned int i=0; i -double Coefficient::value (const Point &/*p*/, - const unsigned int /*component*/) const -{ - return 1; -} @@ -415,8 +413,7 @@ TwoPhaseFlowProblem::TwoPhaseFlowProblem (const unsigned int degree) FE_DGQ(degree), 1), dof_handler (triangulation), n_refinement_steps (5), - time_step (1.0/std::pow(2.0, double(n_refinement_steps))/6), - epsilon(0.05), + time_step (10.0/std::pow(2.0, double(n_refinement_steps))/6), vis (0.2) {} @@ -436,11 +433,11 @@ void TwoPhaseFlowProblem::make_grid_and_dofs () GridGenerator::hyper_cube (triangulation, 0, 1); for (unsigned int f=0; f::faces_per_cell; ++f) - { if (triangulation.begin()->face(f)->center()[0] == 0) + { if (triangulation.begin()->face(f)->center()[0] == 0) triangulation.begin()->face(f)->set_boundary_indicator (1); - if (triangulation.begin()->face(f)->center()[0] == 1) - triangulation.begin()->face(f)->set_boundary_indicator (2); - } + if (triangulation.begin()->face(f)->center()[0] == 1) + triangulation.begin()->face(f)->set_boundary_indicator (2); + } triangulation.refine_global (n_refinement_steps); @@ -527,6 +524,7 @@ void TwoPhaseFlowProblem::make_grid_and_dofs () //at that time, we have the new velocity solved // we can use it to assemble Matrixblock(0,2) + const KInverse<2> k_inverse; template void TwoPhaseFlowProblem::assemble_system () @@ -567,14 +565,11 @@ void TwoPhaseFlowProblem::assemble_system () // matrices. const RightHandSide right_hand_side; const PressureBoundaryValues pressure_boundary_values; - const KInverse k_inverse; - const Coefficient coefficient; std::vector rhs_values (n_q_points); std::vector boundary_values (n_face_q_points); std::vector > k_inverse_values (n_q_points); - std::vector coefficient_values(n_q_points); std::vector > old_solution_values(n_q_points, Vector(dim+2)); std::vector > > old_solution_grads(n_q_points, @@ -604,10 +599,7 @@ void TwoPhaseFlowProblem::assemble_system () right_hand_side.value_list (fe_values.get_quadrature_points(), rhs_values); k_inverse.value_list (fe_values.get_quadrature_points(), - k_inverse_values, - epsilon); - - coefficient.value_list (fe_values.get_quadrature_points(), coefficient_values); + k_inverse_values); for (unsigned int q=0; q::assemble_system () mobility_inverse(old_s,vis) * phi_j_u - div_phi_i_u * phi_j_p - phi_i_p * div_phi_j_u - + coefficient_values[q] * phi_i_s * phi_j_s - ) + + phi_i_s * phi_j_s + ) * fe_values.JxW(q); } @@ -641,7 +633,7 @@ void TwoPhaseFlowProblem::assemble_system () fe_values.JxW(q); } - //here, we compute the boundary values for pressure + //here, we compute the boundary values for pressure for (unsigned int face_no=0; face_no::faces_per_cell; @@ -688,10 +680,10 @@ void TwoPhaseFlowProblem::assemble_system () for (unsigned int i=0; i::assemble_system () // @sect4{The ``InverseMatrix'' class template} - // Everything here is completely same with step-20 + // Everything here is completely same with step-20 @@ -837,7 +829,7 @@ void TwoPhaseFlowProblem::solve () Vector tmp2 (solution.block(2).size()); - // this part is for pressure + // this part is for pressure { m_inverse.vmult (tmp, system_rhs.block(0)); system_matrix.block(1,0).vmult (schur_rhs, tmp); @@ -861,7 +853,8 @@ void TwoPhaseFlowProblem::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; } @@ -881,189 +874,186 @@ void TwoPhaseFlowProblem::solve () m_inverse.vmult (solution.block(0), tmp); } - //This part is for saturation. - // Here are many complicated functions - //which are very similiar with the - //assemble_system() part. - // For DG(0), we have to consider the discontinuty - // of the solution, then as in Introduction, - // compute numerical flux and judge it is in-flow or out-flow. - // After assemble Matrixbloc(2,0) - // , we could compute saturation directly. + //This part is for saturation. + // Here are many complicated functions + //which are very similiar with the + //assemble_system() part. + // For DG(0), we have to consider the discontinuty + // of the solution, then as in Introduction, + // compute numerical flux and judge it is in-flow or out-flow. + // After assemble Matrixbloc(2,0) + // , we could compute saturation directly. { - QGauss quadrature_formula(degree+2); - QGauss face_quadrature_formula(degree+2); - FEValues fe_values (fe, quadrature_formula, - update_values | update_gradients | - update_q_points | update_JxW_values); - FEFaceValues fe_face_values (fe, face_quadrature_formula, - update_values | update_normal_vectors | - update_q_points | update_JxW_values); - FEFaceValues fe_face_values_neighbor (fe, face_quadrature_formula, - update_values); + QGauss quadrature_formula(degree+2); + QGauss face_quadrature_formula(degree+2); + FEValues fe_values (fe, quadrature_formula, + update_values | update_gradients | + update_q_points | update_JxW_values); + FEFaceValues fe_face_values (fe, face_quadrature_formula, + update_values | update_normal_vectors | + update_q_points | update_JxW_values); + FEFaceValues fe_face_values_neighbor (fe, face_quadrature_formula, + update_values); - const unsigned int dofs_per_cell = fe.dofs_per_cell; - const unsigned int n_q_points = quadrature_formula.n_quadrature_points; - const unsigned int n_face_q_points = face_quadrature_formula.n_quadrature_points; + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.n_quadrature_points; + const unsigned int n_face_q_points = face_quadrature_formula.n_quadrature_points; - vfs_out = 0.0; - v_out = 0.0; + vfs_out = 0.0; + v_out = 0.0; - Vector local_rhs (dofs_per_cell); - std::vector > old_solution_values(n_q_points, Vector(dim+2)); - std::vector > old_solution_values_face(n_face_q_points, Vector(dim+2)); - std::vector > old_solution_values_face_neighbor(n_face_q_points, Vector(dim+2)); - std::vector > present_solution_values(n_q_points, Vector(dim+2)); - std::vector > present_solution_values_face(n_face_q_points, Vector(dim+2)); - - std::vector neighbor_saturation (n_face_q_points); - std::vector local_dof_indices (dofs_per_cell); + Vector local_rhs (dofs_per_cell); + std::vector > old_solution_values(n_q_points, Vector(dim+2)); + std::vector > old_solution_values_face(n_face_q_points, Vector(dim+2)); + std::vector > old_solution_values_face_neighbor(n_face_q_points, Vector(dim+2)); + std::vector > present_solution_values(n_q_points, Vector(dim+2)); + std::vector > present_solution_values_face(n_face_q_points, Vector(dim+2)); + + std::vector neighbor_saturation (n_face_q_points); + std::vector local_dof_indices (dofs_per_cell); - typename DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); - for (; cell!=endc; ++cell) - { - local_rhs = 0; - fe_values.reinit (cell); + for (; cell!=endc; ++cell) + { + local_rhs = 0; + fe_values.reinit (cell); - fe_values.get_function_values (old_solution, old_solution_values); - fe_values.get_function_values (solution, present_solution_values); + fe_values.get_function_values (old_solution, old_solution_values); + fe_values.get_function_values (solution, present_solution_values); - for (unsigned int q=0; q present_u; - for (unsigned int d=0; d grad_phi_i_s = extract_grad_s(fe_values, i, q); + for (unsigned int q=0; q present_u; + for (unsigned int d=0; d grad_phi_i_s = extract_grad_s(fe_values, i, q); - local_rhs(i) += ( - time_step *(f_saturation(old_s,vis) * present_u * grad_phi_i_s)+ - old_s * phi_i_s) - * fe_values.JxW(q); - } - //Here is our numerical flux computation - // Finding neighbor as step-12 + local_rhs(i) += ( + time_step *(f_saturation(old_s,vis) * present_u * grad_phi_i_s)+ + old_s * phi_i_s) + * fe_values.JxW(q); + } + //Here is our numerical flux computation + // Finding neighbor as step-12 - for (unsigned int face_no=0; face_no::faces_per_cell;++face_no) - { - fe_face_values.reinit (cell, face_no); - - fe_face_values.get_function_values (old_solution, old_solution_values_face); - fe_face_values.get_function_values (solution, present_solution_values_face); - - if (cell->at_boundary(face_no)) - { - if (cell->face(face_no)->boundary_indicator() == 1) - for (unsigned int q=0;q::active_cell_iterator - neighbor = cell->neighbor(face_no); - const unsigned int - neighbor_face = cell->neighbor_of_neighbor(face_no); - - fe_face_values_neighbor.reinit (neighbor, neighbor_face); + for (unsigned int face_no=0; face_no::faces_per_cell;++face_no) + { + fe_face_values.reinit (cell, face_no); + + fe_face_values.get_function_values (old_solution, old_solution_values_face); + fe_face_values.get_function_values (solution, present_solution_values_face); + + if (cell->at_boundary(face_no)) + { + if (cell->face(face_no)->boundary_indicator() == 1) + for (unsigned int q=0;q::active_cell_iterator + neighbor = cell->neighbor(face_no); + const unsigned int + neighbor_face = cell->neighbor_of_neighbor(face_no); + + fe_face_values_neighbor.reinit (neighbor, neighbor_face); - fe_face_values_neighbor.get_function_values (old_solution, - old_solution_values_face_neighbor); + fe_face_values_neighbor.get_function_values (old_solution, + old_solution_values_face_neighbor); - for (unsigned int q=0;qat_boundary(face_no)) - { - if (cell->face(face_no)->boundary_indicator() ==2 ) - {for (unsigned int q=0;q present_u_face; - for (unsigned int d=0; d= 0); + if (cell->at_boundary(face_no)) + { + if (cell->face(face_no)->boundary_indicator() ==2 ) + {for (unsigned int q=0;q present_u_face; + for (unsigned int d=0; d= 0); - if (is_outflow_q_point == true) - { - for (unsigned int i=0; iget_dof_indices (local_dof_indices); - for (unsigned int i=0; iget_dof_indices (local_dof_indices); + for (unsigned int i=0; i cg (solver_control); cg.solve (system_matrix.block(2,2), solution.block(2), system_rhs.block(2), - PreconditionIdentity()); + 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; - - - - + old_solution = solution; } // @sect4{TwoPhaseFlow::compute_errors} @@ -1092,9 +1082,9 @@ Vector TwoPhaseFlowProblem::evaluate_solution (const Point &point) const { static const MappingQ1 mapping; - // first find the cell in which this point - // is, initialize a quadrature rule with - // it, and then a FEValues object + // first find the cell in which this point + // is, initialize a quadrature rule with + // it, and then a FEValues object const typename DoFHandler::active_cell_iterator cell = GridTools::find_active_cell_around_point (dof_handler, point); @@ -1114,15 +1104,15 @@ TwoPhaseFlowProblem::evaluate_solution (const Point &point) const return u_value[0]; } - //{TwoPhaseFlowProblem::compute_errors} + //{TwoPhaseFlowProblem::compute_errors} - // The compute_errors function is to compute - // error on some euqally spaced fixed points - // use evaluation function to interpret - // solution value at the point - // then output those fixed points' value - // For each mesh, we can compare the output - // to estimate errors. + // The compute_errors function is to compute + // error on some euqally spaced fixed points + // use evaluation function to interpret + // solution value at the point + // then output those fixed points' value + // For each mesh, we can compare the output + // to estimate errors. template void TwoPhaseFlowProblem::compute_errors () const @@ -1194,6 +1184,57 @@ void TwoPhaseFlowProblem::output_results } + +template +void +TwoPhaseFlowProblem::project_back_saturation () +{ + for (unsigned int i=0; i 1) + solution.block(dim)(i) = 1; +} + + + +template +double +TwoPhaseFlowProblem::get_maximal_velocity () const +{ + QGauss quadrature_formula(degree+2); + const unsigned int n_q_points + = quadrature_formula.n_quadrature_points; + + FEValues fe_values (fe, quadrature_formula, + update_values); + std::vector > old_solution_values(n_q_points, Vector(dim+2)); + double max_velocity = 0; + + typename DoFHandler::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 velocity; + for (unsigned int i=0; i::output_results template void TwoPhaseFlowProblem::run () { - std::cout<<"Solving problem in " <::run () unsigned int timestep_number = 1; - for ( double time = time_step; time <=2; time+=time_step, timestep_number++) - { - std::cout<< "Timestep_number = "<< timestep_number<::iterator