From: bangerth Date: Fri, 20 Oct 2006 19:52:38 +0000 (+0000) Subject: A few more things X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=314bdba13cf4ccbc50a9318b317614edeb8d3907;p=dealii-svn.git A few more things git-svn-id: https://svn.dealii.org/trunk@14036 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-21/doc/intro.dox b/deal.II/examples/step-21/doc/intro.dox index 12147df1c7..41b27efc19 100644 --- a/deal.II/examples/step-21/doc/intro.dox +++ b/deal.II/examples/step-21/doc/intro.dox @@ -168,11 +168,10 @@ depends on the previously known $S^n$ and the just computed $\mathbf{u}^{n+1}$. We can then state the problem in weak form as follows, by multiplying each equation with test functions $\mathbf v$, $q$, and $\sigma$ and integrating -terms on each cell $K$ by parts: +terms by parts: \f{eqnarray*} - \sum_K - \left((\mathbf{K}\lambda(S^n))^{-1} \mathbf{u}^{n+1},\mathbf v\right)_K - - (p^{n+1}, \mathbf v)_K &=& + \left((\mathbf{K}\lambda(S^n))^{-1} \mathbf{u}^{n+1},\mathbf v\right)_\Omega - + (p^{n+1}, \nabla\cdot\mathbf v)_\Omega &=& - (p^{n+1}, \mathbf v)_{\partial\Omega} \\ (\nabla \cdot\mathbf{u}^{n+1}, q)_\Omega &=& 0 @@ -247,60 +246,93 @@ fluxes can also be found in step-12. The linear solvers used in this program are a straightforward extension of the ones used in step-20. Essentially, we simply have to extend everything from -two to three solution components. - - $(1)$ We has the three blocks vector $(u,p,S)$ , in which -all the functions are dependent on time. i.e. At each time step we -need project the $solution$ into $old-solution$, using -$old-solution$ to get a new $solution$. -Keep doing this until the last time step; - -At time $t=t^{n+1}$ , suppose $old-solution=(u^{n},p^{n},S^{n})$ is -known, in $assemble-system()$ part, we assemble system matrix as: +two to three solution components. If we use the discrete spaces +mentioned above and put shape functions into the bilinear forms, we +arrive at the following linear system to be solved for time step $n+1$: \f[ -\begin{array}{cccccccccccc} -\lceil &M(S^{n}) &B^{T}& 0 &\rceil & \lceil& \mathbf{u}^{n+1}&\rceil& &\lceil& 0 &\rceil\\ -| &B& 0 & 0 & | &| & p^{n+1} &| &=&| & q &|\\ -\lfloor&\triangle t \nabla F(S^n)& 0& I & \rfloor & \lfloor -&S^{n+1} & \rfloor & & \lfloor& S^{n}& \rfloor +\left( +\begin{array}{ccc} +M^u(S^{n}) & B^{T}& 0\\ +B & 0 & 0\\ +\triangle t\; H & 0& M^S +\end{array} +\right) +\left( +\begin{array}{c} +\mathbf{U}^{n+1} \\ P^{n+1} \\ S^{n+1} +\end{array} +\right) += +\left( +\begin{array}{c} +0 \\ F_2 \\ F_3 \end{array} +\right) \f] +where the individual matrices and vectors are defined as follows using +shape functions $\mathbf v_i$ (of type Raviart Thomas $RT_k$) for +velocities and $\phi_i$ (of type $DG_k$) for both pressures and saturations: +\f{eqnarray*} +M^u(S^n)_{ij} &=& +\left((\mathbf{K}\lambda(S^n))^{-1} \mathbf{v}_i,\mathbf +v_j\right)_\Omega, +\\ +B_{ij} &=& +(\nabla \cdot \mathbf v_i, \phi_j)_\Omega, +\\ +H_{ij} &=& + \sum_K + \left\{ + \left(F(S^n) \mathbf v_i, \nabla \phi_j)\right)_K + - + \left(F(S^n_+) (\mathbf n \cdot (\mathbf v_i)_+), \phi_j\right)_{\partial K_+} + - + \left(F(S^n_-) (\mathbf n \cdot (\mathbf v_i)_-), \phi_j\right)_{\partial K_-}, + \right\} +\\ +M^S_{ij} &=& +(\phi_i, \phi_j)_\Omega, +\\ +(F_2)_i &=& +(q,\phi_i)_\Omega, +\\ +(F_3)_i &=& +(S^n,\phi_i)_\Omega. +\f} -In $solve()$ part, we solve the first two equations independent of -the third equation, since $M( S^n)$ is already known. As in step-20, -using vector base functions, Schur complement with a -preconditioner and CG method, we get $u^{n+1}$and $p^{n+1}$. +Note the following complication, however: Since the matrix $H_{ij}$ +depends on $\mathbf u^{n+1}$ implicitly (the velocities are needed to +determine which parts of the boundaries $\partial K$ of cells are +influx or outflux parts), we can only assemble this matrix after we +have solved for the velocities. -Then, with the above $u^{n+1}$ and $p^{n+1}$, we could compute -$S^{n+1}$ by : -@f{eqnarray*} -\sum_{\kappa}\int_{\kappa}S^{n+1} \phi dx&&= -\sum_{\kappa}\{\int_{\kappa} S^{n} \phi dx+\Delta t\int_{\kappa} -F(S^{n}) \mathbf{u^{n+1}}\cdot \nabla \phi dx\\ -\nonumber && -\Delta t \int_{\partial \kappa_{-}}F(S^{n,-}) -\mathbf{u}^{n+1,-}\cdot \mathbf{n} \phi dx -\Delta t \int_{\partial -\kappa_{+}}F(S^{n}) \mathbf{u}^{n+1}\cdot \mathbf{n} \phi dx\} -@f} -Now, project solution $(u^{n+1},p^{n+1},S^{n+1})$ into -$old-solution$, do the above process for next time step. +The solution scheme then involves the following steps: +
    +
  1. Solve for the pressure $p^{n+1}$ using the Schur complement + technique introduced in step-20. + +
  2. Solve for the velocity $\mathbf u^{n+1}$ as also discussed in + step-20. + +
  3. Compute the term $F_3-\triangle t\; H \mathbf u^{n+1}$, using + the just computed velocities. + +
  4. Solve for the saturation $S^{n+1}$. +
- -$(2)$ The numerical flux term is related with neighbor cells. In our -implementation $solve( )$, we do the following on each cell: +In this scheme, we never actually build the matrix $H$, but rather +generate the right hand side of the third equation once we are ready +to do so. -For each face, compute the flux $\mathbf{u}\cdot F(S)$, the flux is -negative means it is an in-flow face. Then if this in-flow face is -on the boundary $\Gamma_{1}$:$F(S^{-})=F(1)$; If the in-flow -face is not on boundary, $F(S^{-})=F(S|_{neighbor})$. -Flux is positive means it is an out-flow face, we just use $ F(S)$ on -current cell. +In the program, we use a variable solution to store the +solution of the present time step. At the end of each step, we copy +its content, i.e. all three of its block components, into the variable +old_solution for use in the next time step. -All the other functions are commented in code, please see next part -- the commented program.

Test Case

- For simplicity, in our project we will assume no +For simplicity, in our project we will assume no source $q=0$ and the heterogeneous porous medium is isotropic $\mathbf{K}(x,y) = k(x,y) \mathbf{I}$. diff --git a/deal.II/examples/step-21/step-21.cc b/deal.II/examples/step-21/step-21.cc index 4304039460..6c0ab4bcd1 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 @@ -81,7 +81,6 @@ class TwoPhaseFlowProblem const unsigned int n_refinement_steps; double time_step; - double epsilon; double vis; double vfs_out; double v_out; @@ -96,7 +95,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. @@ -196,11 +195,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 @@ -256,68 +255,69 @@ 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())); - + for (unsigned int p=0; p -double Coefficient::value (const Point &/*p*/, - const unsigned int /*component*/) const -{ - return 1; -} @@ -408,8 +408,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) {} @@ -429,11 +428,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); @@ -561,13 +560,11 @@ void TwoPhaseFlowProblem::assemble_system () 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, @@ -597,10 +594,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); } @@ -634,7 +628,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; @@ -681,10 +675,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 @@ -830,7 +824,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); @@ -874,176 +868,176 @@ 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() @@ -1052,7 +1046,7 @@ void TwoPhaseFlowProblem::solve () } - old_solution = solution; + old_solution = solution; @@ -1085,9 +1079,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); @@ -1107,15 +1101,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 @@ -1185,7 +1179,7 @@ void TwoPhaseFlowProblem::output_results std::ofstream output (filename.str().c_str()); data_out.write_gnuplot (output); - //data_out.write_vtk (output); + //data_out.write_vtk (output); } @@ -1214,16 +1208,16 @@ void TwoPhaseFlowProblem::run () unsigned int timestep_number = 1; for ( double time = time_step; time <=1; time+=time_step, timestep_number++) - { - std::cout<< "Timestep_number = "<< timestep_number<::iterator