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
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:
+<ol>
+ <li>Solve for the pressure $p^{n+1}$ using the Schur complement
+ technique introduced in step-20.
+
+ <li>Solve for the velocity $\mathbf u^{n+1}$ as also discussed in
+ step-20.
+
+ <li>Compute the term $F_3-\triangle t\; H \mathbf u^{n+1}$, using
+ the just computed velocities.
+
+ <li>Solve for the saturation $S^{n+1}$.
+</ol>
-
-$(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 <code>solution</code> 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
+<code>old_solution</code> for use in the next time step.
-All the other functions are commented in code, please see next part
-- the commented program.
<h2>Test Case</h2>
- 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}$.
- // @sect3{Include files}
+ // @sect3{Include files}
// This program is an daptation of step-20
// and includes some technique of DG method from step-12
#include <dofs/dof_tools.h>
#include <dofs/dof_constraints.h>
#include <fe/fe_q.h>
- //The Discontinuous Galerkin finite element is declared:
+ //The Discontinuous Galerkin finite element is declared:
#include <fe/fe_dgq.h>
#include <fe/fe_system.h>
#include <fstream>
#include <iostream>
#include <sstream>
- // The Raviart-Thomas finite element is declared:
+ // The Raviart-Thomas finite element is declared:
#include <fe/fe_raviart_thomas.h>
// In this program, we use a tensorial
const unsigned int n_refinement_steps;
double time_step;
- double epsilon;
double vis;
double vfs_out;
double v_out;
//{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.
template <int dim>
double InitialValues<dim>::value (const Point<dim> &p,
- const unsigned int component) const
+ const unsigned int component) const
{
- if(component<dim+1)
+ if(component<dim+1)
return 0;
- else
+ else
{
if(p[0]==0)return 1;
else return 0;
// ``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
// of the function in the second
// argument, a list of tensors:
template <int dim>
-class KInverse //: public TensorFunction<2,dim>
+class KInverse : public TensorFunction<2,dim>
{
public:
- KInverse ()
- {}
+ KInverse ();
- /*virtual*/ void value_list (const std::vector<Point<dim> > &points,
- std::vector<Tensor<2,dim> > &values,
- const double epsilon) const;
+ virtual void value_list (const std::vector<Point<dim> > &points,
+ std::vector<Tensor<2,dim> > &values) const;
+
+ private:
+ std::vector<Point<dim> > centers;
};
template <int dim>
-class Coefficient : public Function<dim>
+KInverse<dim>::KInverse ()
{
- public:
- Coefficient () : Function<dim>(1) {};
-
- virtual double value (const Point<dim> &p,
- const unsigned int component = 0 ) const;
-};
+ const unsigned int N = 40;
+ centers.resize (N);
+ for (unsigned int i=0; i<N; ++i)
+ for (unsigned int d=0; d<dim; ++d)
+ centers[i][d] = 2.*rand()/RAND_MAX-1;
+}
+
template <int dim>
void
KInverse<dim>::value_list (const std::vector<Point<dim> > &points,
- std::vector<Tensor<2,dim> > &values,
- const double epsilon) const
+ std::vector<Tensor<2,dim> > &values) const
{
Assert (points.size() == values.size(),
ExcDimensionMismatch (points.size(), values.size()));
-
+
for (unsigned int p=0; p<points.size(); ++p)
{
values[p].clear ();
- const double permeability = 1.0/(2+1.99*sin(2*3.1415926*(2*points[p][0]-1*points[p][1])/epsilon));
+ double permeability = 0;
+ for (unsigned int i=0; i<centers.size(); ++i)
+ permeability += std::exp(-(points[p]-centers[i]).square()
+ / (0.1 * 0.1));
+
+ const double normalized_permeability
+ = std::max(permeability, 0.005);
for (unsigned int d=0; d<dim; ++d)
- values[p][d][d]=1.0/permeability;
-
+ values[p][d][d] = 10./normalized_permeability;
}
}
+
double mobility_inverse (const double S, const double vis)
{
- return 1.0 /(1.0/vis * S * S + (1-S) * (1-S));
+ 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));
+ return S*S /( S * S +vis * (1-S) * (1-S));
}
-template <int dim>
-double Coefficient<dim>::value (const Point<dim> &/*p*/,
- const unsigned int /*component*/) const
-{
- return 1;
-}
FE_DGQ<dim>(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)
{}
GridGenerator::hyper_cube (triangulation, 0, 1);
for (unsigned int f=0; f<GeometryInfo<dim>::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);
const RightHandSide<dim> right_hand_side;
const PressureBoundaryValues<dim> pressure_boundary_values;
const KInverse<dim> k_inverse;
- const Coefficient<dim> coefficient;
std::vector<double> rhs_values (n_q_points);
std::vector<double> boundary_values (n_face_q_points);
std::vector<Tensor<2,dim> > k_inverse_values (n_q_points);
- std::vector<double> coefficient_values(n_q_points);
std::vector<Vector<double> > old_solution_values(n_q_points, Vector<double>(dim+2));
std::vector<std::vector<Tensor<1,dim> > > old_solution_grads(n_q_points,
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<n_q_points; ++q)
for (unsigned int i=0; i<dofs_per_cell; ++i)
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);
}
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<GeometryInfo<dim>::faces_per_cell;
for (unsigned int i=0; i<dofs_per_cell; ++i)
for (unsigned int j=0; j<dofs_per_cell; ++j)
- { system_matrix.add (local_dof_indices[i],
- local_dof_indices[j],
- local_matrix(i,j));
- }
+ { system_matrix.add (local_dof_indices[i],
+ local_dof_indices[j],
+ local_matrix(i,j));
+ }
for (unsigned int i=0; i<dofs_per_cell; ++i)
system_rhs(local_dof_indices[i]) += local_rhs(i);
// @sect4{The ``InverseMatrix'' class template}
- // Everything here is completely same with step-20
+ // Everything here is completely same with step-20
Vector<double> 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);
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<dim> quadrature_formula(degree+2);
- QGauss<dim-1> face_quadrature_formula(degree+2);
- FEValues<dim> fe_values (fe, quadrature_formula,
- update_values | update_gradients |
- update_q_points | update_JxW_values);
- FEFaceValues<dim> fe_face_values (fe, face_quadrature_formula,
- update_values | update_normal_vectors |
- update_q_points | update_JxW_values);
- FEFaceValues<dim> fe_face_values_neighbor (fe, face_quadrature_formula,
- update_values);
+ QGauss<dim> quadrature_formula(degree+2);
+ QGauss<dim-1> face_quadrature_formula(degree+2);
+ FEValues<dim> fe_values (fe, quadrature_formula,
+ update_values | update_gradients |
+ update_q_points | update_JxW_values);
+ FEFaceValues<dim> fe_face_values (fe, face_quadrature_formula,
+ update_values | update_normal_vectors |
+ update_q_points | update_JxW_values);
+ FEFaceValues<dim> 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<double> local_rhs (dofs_per_cell);
- std::vector<Vector<double> > old_solution_values(n_q_points, Vector<double>(dim+2));
- std::vector<Vector<double> > old_solution_values_face(n_face_q_points, Vector<double>(dim+2));
- std::vector<Vector<double> > old_solution_values_face_neighbor(n_face_q_points, Vector<double>(dim+2));
- std::vector<Vector<double> > present_solution_values(n_q_points, Vector<double>(dim+2));
- std::vector<Vector<double> > present_solution_values_face(n_face_q_points, Vector<double>(dim+2));
-
- std::vector<double> neighbor_saturation (n_face_q_points);
- std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+ Vector<double> local_rhs (dofs_per_cell);
+ std::vector<Vector<double> > old_solution_values(n_q_points, Vector<double>(dim+2));
+ std::vector<Vector<double> > old_solution_values_face(n_face_q_points, Vector<double>(dim+2));
+ std::vector<Vector<double> > old_solution_values_face_neighbor(n_face_q_points, Vector<double>(dim+2));
+ std::vector<Vector<double> > present_solution_values(n_q_points, Vector<double>(dim+2));
+ std::vector<Vector<double> > present_solution_values_face(n_face_q_points, Vector<double>(dim+2));
+
+ std::vector<double> neighbor_saturation (n_face_q_points);
+ std::vector<unsigned int> local_dof_indices (dofs_per_cell);
- typename DoFHandler<dim>::active_cell_iterator
- cell = dof_handler.begin_active(),
- endc = dof_handler.end();
+ typename DoFHandler<dim>::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<n_q_points; ++q)
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- {
- const double old_s = old_solution_values[q](dim+1);
- Tensor<1,dim> present_u;
- for (unsigned int d=0; d<dim; ++d)
- present_u[d] = present_solution_values[q](d);
-
- const double phi_i_s = extract_s(fe_values, i, q);
- const Tensor<1,dim> grad_phi_i_s = extract_grad_s(fe_values, i, q);
+ for (unsigned int q=0; q<n_q_points; ++q)
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ const double old_s = old_solution_values[q](dim+1);
+ Tensor<1,dim> present_u;
+ for (unsigned int d=0; d<dim; ++d)
+ present_u[d] = present_solution_values[q](d);
+
+ const double phi_i_s = extract_s(fe_values, i, q);
+ const Tensor<1,dim> 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<GeometryInfo<dim>::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<n_face_q_points;++q)
- neighbor_saturation[q] = 1;
- else
- for (unsigned int q=0;q<n_face_q_points;++q)
- neighbor_saturation[q] = 0;
- }
- else
- // there is a neighbor behind this face
- {
- const typename DoFHandler<dim>::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<GeometryInfo<dim>::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<n_face_q_points;++q)
+ neighbor_saturation[q] = 1;
+ else
+ for (unsigned int q=0;q<n_face_q_points;++q)
+ neighbor_saturation[q] = 0;
+ }
+ else
+ // there is a neighbor behind this face
+ {
+ const typename DoFHandler<dim>::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;q<n_face_q_points;++q)
- neighbor_saturation[q] = old_solution_values_face_neighbor[q](dim+1);
- }
+ for (unsigned int q=0;q<n_face_q_points;++q)
+ neighbor_saturation[q] = old_solution_values_face_neighbor[q](dim+1);
+ }
- if (cell->at_boundary(face_no))
- {
- if (cell->face(face_no)->boundary_indicator() ==2 )
- {for (unsigned int q=0;q<n_face_q_points;++q)
- {
- vfs_out += present_solution_values_face[q](0)
- *f_saturation(present_solution_values_face[q](dim+1),vis)
- *fe_face_values.JxW(q);
- v_out += present_solution_values_face[q](0)
- *fe_face_values.JxW(q);
- }
- }
- }
- for (unsigned int q=0;q<n_face_q_points;++q)
- {
- Tensor<1,dim> present_u_face;
- for (unsigned int d=0; d<dim; ++d)
- { present_u_face[d] = present_solution_values_face[q](d);
- }
- const double normal_flux = present_u_face *
- fe_face_values.normal_vector(q);
-
- const bool is_outflow_q_point = (normal_flux >= 0);
+ if (cell->at_boundary(face_no))
+ {
+ if (cell->face(face_no)->boundary_indicator() ==2 )
+ {for (unsigned int q=0;q<n_face_q_points;++q)
+ {
+ vfs_out += present_solution_values_face[q](0)
+ *f_saturation(present_solution_values_face[q](dim+1),vis)
+ *fe_face_values.JxW(q);
+ v_out += present_solution_values_face[q](0)
+ *fe_face_values.JxW(q);
+ }
+ }
+ }
+ for (unsigned int q=0;q<n_face_q_points;++q)
+ {
+ Tensor<1,dim> present_u_face;
+ for (unsigned int d=0; d<dim; ++d)
+ { present_u_face[d] = present_solution_values_face[q](d);
+ }
+ const double normal_flux = present_u_face *
+ fe_face_values.normal_vector(q);
+
+ const bool is_outflow_q_point = (normal_flux >= 0);
- if (is_outflow_q_point == true)
- {
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- {
- const double outflow = -time_step * normal_flux
- * f_saturation(old_solution_values_face[q](dim+1),vis)
- * extract_s(fe_face_values,i,q)
- * fe_face_values.JxW(q);
- local_rhs(i) += outflow;
- }
- }
+ if (is_outflow_q_point == true)
+ {
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ const double outflow = -time_step * normal_flux
+ * f_saturation(old_solution_values_face[q](dim+1),vis)
+ * extract_s(fe_face_values,i,q)
+ * fe_face_values.JxW(q);
+ local_rhs(i) += outflow;
+ }
+ }
- else
- {
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- {
- const double inflow = -time_step * normal_flux
- * f_saturation( neighbor_saturation[q],vis)
- * extract_s(fe_face_values,i,q)
- * fe_face_values.JxW(q);
- local_rhs(i) += inflow;
- }
+ else
+ {
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ const double inflow = -time_step * normal_flux
+ * f_saturation( neighbor_saturation[q],vis)
+ * extract_s(fe_face_values,i,q)
+ * fe_face_values.JxW(q);
+ local_rhs(i) += inflow;
+ }
- }
+ }
- }
+ }
- }
+ }
- cell->get_dof_indices (local_dof_indices);
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- {
- system_rhs(local_dof_indices[i]) += local_rhs(i);
- }
+ cell->get_dof_indices (local_dof_indices);
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ system_rhs(local_dof_indices[i]) += local_rhs(i);
+ }
- }
+ }
SolverControl solver_control (system_matrix.block(2,2).m(),
1e-12*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());
+ PreconditionIdentity());
std::cout << solver_control.last_step()
}
- old_solution = solution;
+ old_solution = solution;
TwoPhaseFlowProblem<dim>::evaluate_solution (const Point<dim> &point) const
{
static const MappingQ1<dim> 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<dim>::active_cell_iterator
cell = GridTools::find_active_cell_around_point (dof_handler, point);
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 <int dim>
void TwoPhaseFlowProblem<dim>::compute_errors () const
std::ofstream output (filename.str().c_str());
data_out.write_gnuplot (output);
- //data_out.write_vtk (output);
+ //data_out.write_vtk (output);
}
unsigned int timestep_number = 1;
for ( double time = time_step; time <=1; time+=time_step, timestep_number++)
- {
- std::cout<< "Timestep_number = "<< timestep_number<<std::endl;
- assemble_system ();
- solve ();
- 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<< "Timestep_number = "<< timestep_number<<std::endl;
+ assemble_system ();
+ solve ();
+ 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::ofstream production_history ("production_history");
std::list<double>::iterator