From: bangerth Date: Tue, 23 Nov 2010 04:34:25 +0000 (+0000) Subject: Add the actual file. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=6892c25e332c3b77880b3c33c536e6f8c8f34c79;p=dealii-svn.git Add the actual file. git-svn-id: https://svn.dealii.org/trunk@22845 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-43/step-43.cc b/deal.II/examples/step-43/step-43.cc index 3a2cfde492..de380e660a 100644 --- a/deal.II/examples/step-43/step-43.cc +++ b/deal.II/examples/step-43/step-43.cc @@ -11,5 +11,1896 @@ /* to the file deal.II/doc/license.html for the text and */ /* further information on this license. */ -int main () + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include + +using namespace dealii; + +namespace LinearSolvers +{ + template + class InverseMatrix : public Subscriptor + { + public: + InverseMatrix (const Matrix &m, + const Preconditioner &preconditioner); + + + template + void vmult (VectorType &dst, + const VectorType &src) const; + + private: + const SmartPointer matrix; + const Preconditioner &preconditioner; + }; + + + template + InverseMatrix:: + InverseMatrix (const Matrix &m, + const Preconditioner &preconditioner) + : + matrix (&m), + preconditioner (preconditioner) + {} + + + + template + template + void + InverseMatrix:: + vmult (VectorType &dst, + const VectorType &src) const + { + SolverControl solver_control (src.size(), 1e-7*src.l2_norm()); + SolverCG cg (solver_control); + + dst = 0; + + try + { + cg.solve (*matrix, dst, src, preconditioner); + } + catch (std::exception &e) + { + Assert (false, ExcMessage(e.what())); + } + } + + template + class BlockSchurPreconditioner : public Subscriptor + { + public: + BlockSchurPreconditioner ( + const TrilinosWrappers::BlockSparseMatrix &S, + const InverseMatrix &Mpinv, + const PreconditionerA &Apreconditioner); + + void vmult (TrilinosWrappers::BlockVector &dst, + const TrilinosWrappers::BlockVector &src) const; + + private: + const SmartPointer darcy_matrix; + const SmartPointer > m_inverse; + const PreconditionerA &a_preconditioner; + + mutable TrilinosWrappers::Vector tmp; + }; + + + + template + BlockSchurPreconditioner:: + BlockSchurPreconditioner(const TrilinosWrappers::BlockSparseMatrix &S, + const InverseMatrix &Mpinv, + const PreconditionerA &Apreconditioner) + : + darcy_matrix (&S), + m_inverse (&Mpinv), + a_preconditioner (Apreconditioner), + tmp (darcy_matrix->block(1,1).m()) + {} + + + template + void BlockSchurPreconditioner::vmult ( + TrilinosWrappers::BlockVector &dst, + const TrilinosWrappers::BlockVector &src) const + { + a_preconditioner.vmult (dst.block(0), src.block(0)); + darcy_matrix->block(1,0).residual(tmp, dst.block(0), src.block(1)); + tmp *= -1; + m_inverse->vmult (dst.block(1), tmp); + } +} + + +template +class TwoPhaseFlowProblem +{ + public: + TwoPhaseFlowProblem (const unsigned int degree); + void run (); + + private: + void setup_dofs (); + void assemble_darcy_preconditioner (); + void build_darcy_preconditioner (); + void assemble_darcy_system (); + void assemble_saturation_system (); + void assemble_saturation_matrix (); + void assemble_saturation_rhs (); + void assemble_saturation_rhs_cell_term (const FEValues &saturation_fe_values, + const FEValues &darcy_fe_values, + const std::vector &local_dof_indices, + const double global_u_infty, + const double global_S_variation, + const double global_Omega_diameter); + void assemble_saturation_rhs_boundary_term (const FEFaceValues &saturation_fe_face_values, + const FEFaceValues &darcy_fe_face_values, + const std::vector &local_dof_indices); + double get_maximal_velocity () const; + std::pair get_extrapolated_saturation_range () const; + void solve (); + bool determine_whether_to_solve_pressure_velocity_part () const; + void compute_refinement_indicators (Vector &indicator) const; + void refine_grid (const Vector &indicator); + void project_back_saturation (); + void output_results () const; + + static + double + compute_viscosity(const std::vector &old_saturation, + const std::vector &old_old_saturation, + const std::vector > &old_saturation_grads, + const std::vector > &old_old_saturation_grads, + const std::vector > &present_darcy_values, + const double global_u_infty, + const double global_S_variation, + const double global_Omega_diameter, + const double cell_diameter, + const double old_time_step, + const double viscosity); + + + const unsigned int degree; + + Triangulation triangulation; + + const unsigned int darcy_degree; + FESystem darcy_fe; + DoFHandler darcy_dof_handler; + ConstraintMatrix darcy_constraints; + + ConstraintMatrix darcy_preconditioner_constraints; + + TrilinosWrappers::BlockSparseMatrix darcy_matrix; + TrilinosWrappers::BlockSparseMatrix darcy_preconditioner_matrix; + + TrilinosWrappers::BlockVector darcy_solution; + TrilinosWrappers::BlockVector darcy_rhs; + + TrilinosWrappers::BlockVector nth_darcy_solution_after_solving_pressure_part; + TrilinosWrappers::BlockVector n_minus_oneth_darcy_solution_after_solving_pressure_part; + + const unsigned int saturation_degree; + FE_Q saturation_fe; + DoFHandler saturation_dof_handler; + ConstraintMatrix saturation_constraints; + + TrilinosWrappers::SparseMatrix saturation_matrix; + + TrilinosWrappers::Vector predictor_saturation_solution; + TrilinosWrappers::Vector saturation_solution; + TrilinosWrappers::Vector old_saturation_solution; + TrilinosWrappers::Vector old_old_saturation_solution; + TrilinosWrappers::Vector saturation_rhs; + + TrilinosWrappers::Vector nth_saturation_solution_after_solving_pressure_part; + + const unsigned int n_refinement_steps; + bool solve_pressure_velocity_part; + bool previous_solve_pressure_velocity_part; + + const double saturation_level; + const double saturation_value; + + double n_minus_oneth_time_step; + double cumulative_nth_time_step; + + double time_step; + double old_time_step; + unsigned int timestep_number; + double viscosity; + + std_cxx1x::shared_ptr Amg_preconditioner; + std_cxx1x::shared_ptr Mp_preconditioner; + + bool rebuild_saturation_matrix; +}; + + +template +class PressureRightHandSide : public Function +{ + public: + PressureRightHandSide () : Function(1) {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; +}; + + + +template +double +PressureRightHandSide::value (const Point &/*p*/, + const unsigned int /*component*/) const +{ + return 0; +} + + +template +class PressureBoundaryValues : public Function +{ + public: + PressureBoundaryValues () : Function(1) {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; +}; + + +template +double +PressureBoundaryValues::value (const Point &p, + const unsigned int /*component*/) const +{ + return 1-p[0]; +} + + +template +class SaturationBoundaryValues : public Function +{ + public: + SaturationBoundaryValues () : Function(1) {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; +}; + + + +template +double +SaturationBoundaryValues::value (const Point &p, + const unsigned int /*component*/) const +{ + if (p[0] == 0) + return 1; + else + return 0; +} + + +template +class SaturationInitialValues : public Function +{ + public: + SaturationInitialValues () : Function(1) {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; + + virtual void vector_value (const Point &p, + Vector &value) const; + +}; + + +template +double +SaturationInitialValues::value (const Point &/*p*/, + const unsigned int /*component*/) const +{ + return 0; +} + + +template +void +SaturationInitialValues::vector_value (const Point &p, + Vector &values) const +{ + for (unsigned int c=0; cn_components; ++c) + values(c) = SaturationInitialValues::value (p,c); +} + + +namespace SingleCurvingCrack +{ + template + class KInverse : public TensorFunction<2,dim> + { + public: + KInverse () + : + TensorFunction<2,dim> () + {} + + virtual void value_list (const std::vector > &points, + std::vector > &values) const; + }; + + + template + void + KInverse::value_list (const std::vector > &points, + std::vector > &values) const + { + Assert (points.size() == values.size(), + ExcDimensionMismatch (points.size(), values.size())); + + for (unsigned int p=0; p + class KInverse : public TensorFunction<2,dim> + { + public: + KInverse () + : + TensorFunction<2,dim> () + {} + + virtual void value_list (const std::vector > &points, + std::vector > &values) const; + + private: + static std::vector > centers; + + static std::vector > get_centers (); + }; + + + + template + std::vector > + KInverse::centers = KInverse::get_centers(); + + + template + std::vector > + KInverse::get_centers () + { + const unsigned int N = (dim == 2 ? + 40 : + (dim == 3 ? + 100 : + throw ExcNotImplemented())); + + std::vector > centers_list (N); + for (unsigned int i=0; i(rand())/RAND_MAX; + + return centers_list; + } + + + + template + void + KInverse::value_list (const std::vector > &points, + std::vector > &values) const + { + Assert (points.size() == values.size(), + ExcDimensionMismatch (points.size(), values.size())); + + for (unsigned int p=0; p +TwoPhaseFlowProblem::TwoPhaseFlowProblem (const unsigned int degree) + : + degree (degree), + darcy_degree (degree), + darcy_fe (FE_Q(darcy_degree+1), dim, + FE_Q(darcy_degree), 1), + darcy_dof_handler (triangulation), + + saturation_degree (degree), + saturation_fe (saturation_degree), + saturation_dof_handler (triangulation), + + n_refinement_steps (4), + solve_pressure_velocity_part (false), + previous_solve_pressure_velocity_part (false), + + saturation_level (2), + saturation_value (0.5), + + time_step (0), + old_time_step (0), + viscosity (0.2), + + rebuild_saturation_matrix (true) {} + + +template +void TwoPhaseFlowProblem::setup_dofs () +{ + std::vector darcy_block_component (dim+1,0); + darcy_block_component[dim] = 1; + { + darcy_dof_handler.distribute_dofs (darcy_fe); + DoFRenumbering::Cuthill_McKee (darcy_dof_handler); + DoFRenumbering::component_wise (darcy_dof_handler, darcy_block_component); + + darcy_constraints.clear (); + DoFTools::make_hanging_node_constraints (darcy_dof_handler, darcy_constraints); + darcy_constraints.close (); + } + { + saturation_dof_handler.distribute_dofs (saturation_fe); + + saturation_constraints.clear (); + DoFTools::make_hanging_node_constraints (saturation_dof_handler, saturation_constraints); + saturation_constraints.close (); + } + { + darcy_preconditioner_constraints.clear (); + + std::vector component_mask (dim+1, false); + component_mask[dim] = true; + + + DoFTools::make_hanging_node_constraints (darcy_dof_handler, darcy_preconditioner_constraints); + DoFTools::make_zero_boundary_constraints (darcy_dof_handler, darcy_preconditioner_constraints, component_mask); + + darcy_preconditioner_constraints.close (); + } + + + std::vector darcy_dofs_per_block (2); + DoFTools::count_dofs_per_block (darcy_dof_handler, darcy_dofs_per_block, darcy_block_component); + const unsigned int n_u = darcy_dofs_per_block[0], + n_p = darcy_dofs_per_block[1], + n_s = saturation_dof_handler.n_dofs(); + + std::cout << "Number of active cells: " + << triangulation.n_active_cells() + << " (on " + << triangulation.n_levels() + << " levels)" + << std::endl + << "Number of degrees of freedom: " + << n_u + n_p + n_s + << " (" << n_u << '+' << n_p << '+'<< n_s <<')' + << std::endl + << std::endl; + + { + darcy_matrix.clear (); + + BlockCompressedSimpleSparsityPattern csp (2,2); + + csp.block(0,0).reinit (n_u, n_u); + csp.block(0,1).reinit (n_u, n_p); + csp.block(1,0).reinit (n_p, n_u); + csp.block(1,1).reinit (n_p, n_p); + + csp.collect_sizes (); + + Table<2,DoFTools::Coupling> coupling (dim+1, dim+1); + + for (unsigned int c=0; c coupling (dim+1, dim+1); + for (unsigned int c=0; c +void +TwoPhaseFlowProblem::assemble_darcy_preconditioner () +{ + std::cout << " Rebuilding darcy preconditioner..." << std::endl; + + darcy_preconditioner_matrix = 0; + + const QGauss quadrature_formula(darcy_degree+2); + FEValues darcy_fe_values (darcy_fe, quadrature_formula, + update_JxW_values | + update_values | + update_gradients | + update_quadrature_points); + FEValues saturation_fe_values (saturation_fe, quadrature_formula, + update_values); + + const unsigned int dofs_per_cell = darcy_fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + const RandomMedium::KInverse k_inverse; +// const SingleCurvingCrack::KInverse k_inverse; + + std::vector > k_inverse_values (n_q_points); + Tensor<2,dim> k_value; + + std::vector old_saturation_values (n_q_points); + + FullMatrix local_matrix (dofs_per_cell, dofs_per_cell); + std::vector local_dof_indices (dofs_per_cell); + + std::vector > phi_u (dofs_per_cell); + std::vector > grad_phi_p (dofs_per_cell); + + const FEValuesExtractors::Vector velocities (0); + const FEValuesExtractors::Scalar pressure (dim); + + typename DoFHandler::active_cell_iterator + cell = darcy_dof_handler.begin_active(), + endc = darcy_dof_handler.end(); + typename DoFHandler::active_cell_iterator + saturation_cell = saturation_dof_handler.begin_active(); + + for (; cell!=endc; ++cell, ++saturation_cell) + { + darcy_fe_values.reinit (cell); + saturation_fe_values.reinit (saturation_cell); + + local_matrix = 0; + + saturation_fe_values.get_function_values (old_saturation_solution, old_saturation_values); + + k_inverse.value_list (darcy_fe_values.get_quadrature_points(), + k_inverse_values); + + for (unsigned int q=0; qget_dof_indices (local_dof_indices); + darcy_preconditioner_constraints.distribute_local_to_global (local_matrix, + local_dof_indices, + darcy_preconditioner_matrix); + } +} + +template +void +TwoPhaseFlowProblem::build_darcy_preconditioner () +{ + assemble_darcy_preconditioner (); + + Amg_preconditioner = std_cxx1x::shared_ptr + (new TrilinosWrappers::PreconditionIC()); + Amg_preconditioner->initialize(darcy_preconditioner_matrix.block(0,0)); + + Mp_preconditioner = std_cxx1x::shared_ptr + (new TrilinosWrappers::PreconditionIC()); + Mp_preconditioner->initialize(darcy_preconditioner_matrix.block(1,1)); + +} + +template +void TwoPhaseFlowProblem::assemble_darcy_system () +{ + darcy_matrix = 0; + darcy_rhs = 0; + + QGauss quadrature_formula(darcy_degree+2); + QGauss face_quadrature_formula(darcy_degree+2); + + FEValues darcy_fe_values (darcy_fe, quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + FEValues saturation_fe_values (saturation_fe, quadrature_formula, + update_values); + + FEFaceValues darcy_fe_face_values (darcy_fe, face_quadrature_formula, + update_values | update_normal_vectors | + update_quadrature_points | update_JxW_values); + + const unsigned int dofs_per_cell = darcy_fe.dofs_per_cell; + + const unsigned int n_q_points = quadrature_formula.size(); + const unsigned int n_face_q_points = face_quadrature_formula.size(); + + FullMatrix local_matrix (dofs_per_cell, dofs_per_cell); + Vector local_rhs (dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + const PressureRightHandSide pressure_right_hand_side; + const PressureBoundaryValues pressure_boundary_values; + const RandomMedium::KInverse k_inverse; +// const SingleCurvingCrack::KInverse k_inverse; + + std::vector pressure_rhs_values (n_q_points); + std::vector boundary_values (n_face_q_points); + std::vector > k_inverse_values (n_q_points); + + std::vector old_saturation_values (n_q_points); + + std::vector > phi_u (dofs_per_cell); + std::vector div_phi_u (dofs_per_cell); + std::vector phi_p (dofs_per_cell); + + const FEValuesExtractors::Vector velocities (0); + const FEValuesExtractors::Scalar pressure (dim); + + typename DoFHandler::active_cell_iterator + cell = darcy_dof_handler.begin_active(), + endc = darcy_dof_handler.end(); + typename DoFHandler::active_cell_iterator + saturation_cell = saturation_dof_handler.begin_active(); + + for (; cell!=endc; ++cell, ++saturation_cell) + { + darcy_fe_values.reinit (cell); + saturation_fe_values.reinit (saturation_cell); + + local_matrix = 0; + local_rhs = 0; + + saturation_fe_values.get_function_values (old_saturation_solution, old_saturation_values); + + pressure_right_hand_side.value_list (darcy_fe_values.get_quadrature_points(), + pressure_rhs_values); + k_inverse.value_list (darcy_fe_values.get_quadrature_points(), + k_inverse_values); + + for (unsigned int q=0; q::faces_per_cell; + ++face_no) + if (cell->at_boundary(face_no)) + { + darcy_fe_face_values.reinit (cell, face_no); + + pressure_boundary_values + .value_list (darcy_fe_face_values.get_quadrature_points(), + boundary_values); + + for (unsigned int q=0; q + phi_i_u = darcy_fe_face_values[velocities].value (i, q); + + local_rhs(i) += -(phi_i_u * + darcy_fe_face_values.normal_vector(q) * + boundary_values[q] * + darcy_fe_face_values.JxW(q)); + } + } + + for (unsigned int i=0; iget_dof_indices (local_dof_indices); + + darcy_constraints.distribute_local_to_global (local_matrix, + local_rhs, + local_dof_indices, + darcy_matrix, + darcy_rhs); + + } +} + +template +void TwoPhaseFlowProblem::assemble_saturation_system () +{ + if ( rebuild_saturation_matrix == true ) + { + saturation_matrix = 0; + assemble_saturation_matrix (); + } + + saturation_rhs = 0; + assemble_saturation_rhs (); +} + +template +void TwoPhaseFlowProblem::assemble_saturation_matrix () +{ + QGauss quadrature_formula(saturation_degree+2); + + FEValues saturation_fe_values (saturation_fe, quadrature_formula, + update_values | update_JxW_values); + + const unsigned int dofs_per_cell = saturation_fe.dofs_per_cell; + + const unsigned int n_q_points = quadrature_formula.size(); + + FullMatrix local_matrix (dofs_per_cell, dofs_per_cell); + Vector local_rhs (dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + typename DoFHandler::active_cell_iterator + cell = saturation_dof_handler.begin_active(), + endc = saturation_dof_handler.end(); + for (; cell!=endc; ++cell) + { + saturation_fe_values.reinit (cell); + local_matrix = 0; + local_rhs = 0; + + for (unsigned int q=0; qget_dof_indices (local_dof_indices); + + saturation_constraints.distribute_local_to_global (local_matrix, + local_dof_indices, + saturation_matrix); + + } +} + + +template +void TwoPhaseFlowProblem::assemble_saturation_rhs () +{ + QGauss quadrature_formula(saturation_degree+2); + QGauss face_quadrature_formula(saturation_degree+2); + + FEValues saturation_fe_values (saturation_fe, quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + FEValues darcy_fe_values (darcy_fe, quadrature_formula, + update_values); + FEFaceValues saturation_fe_face_values (saturation_fe, face_quadrature_formula, + update_values | update_normal_vectors | + update_quadrature_points | update_JxW_values); + FEFaceValues darcy_fe_face_values (darcy_fe, face_quadrature_formula, + update_values); + FEFaceValues saturation_fe_face_values_neighbor (saturation_fe, face_quadrature_formula, + update_values); + + const unsigned int dofs_per_cell = saturation_dof_handler.get_fe().dofs_per_cell; + std::vector local_dof_indices (dofs_per_cell); + + const double global_u_infty = get_maximal_velocity (); + const std::pair + global_S_range = get_extrapolated_saturation_range (); + const double global_S_variasion = global_S_range.second - global_S_range.first; + const double global_Omega_diameter = GridTools::diameter (triangulation); + + typename DoFHandler::active_cell_iterator + cell = saturation_dof_handler.begin_active(), + endc = saturation_dof_handler.end(); + typename DoFHandler::active_cell_iterator + darcy_cell = darcy_dof_handler.begin_active(); + for (; cell!=endc; ++cell, ++darcy_cell) + { + saturation_fe_values.reinit (cell); + darcy_fe_values.reinit (darcy_cell); + + cell->get_dof_indices (local_dof_indices); + + assemble_saturation_rhs_cell_term(saturation_fe_values, + darcy_fe_values, + local_dof_indices, + global_u_infty, + global_S_variasion, + global_Omega_diameter); + + for (unsigned int face_no=0; face_no::faces_per_cell; + ++face_no) + { + + if (cell->at_boundary(face_no)) + { + darcy_fe_face_values.reinit (darcy_cell, face_no); + saturation_fe_face_values.reinit (cell, face_no); + assemble_saturation_rhs_boundary_term (saturation_fe_face_values, + darcy_fe_face_values, + local_dof_indices); + } + } + } +} + +template +void +TwoPhaseFlowProblem:: +assemble_saturation_rhs_cell_term (const FEValues &saturation_fe_values, + const FEValues &darcy_fe_values, + const std::vector &local_dof_indices, + const double global_u_infty, + const double global_S_variation, + const double global_Omega_diameter) +{ + const unsigned int dofs_per_cell = saturation_fe_values.dofs_per_cell; + const unsigned int n_q_points = saturation_fe_values.n_quadrature_points; + + Vector local_rhs (dofs_per_cell); + + std::vector old_saturation_solution_values(n_q_points); + std::vector old_old_saturation_solution_values(n_q_points); + std::vector > old_grad_saturation_solution_values(n_q_points); + std::vector > old_old_grad_saturation_solution_values(n_q_points); + std::vector > present_darcy_solution_values(n_q_points, Vector(dim+1)); + + saturation_fe_values.get_function_values (old_saturation_solution, old_saturation_solution_values); + saturation_fe_values.get_function_values (old_old_saturation_solution, old_old_saturation_solution_values); + saturation_fe_values.get_function_grads (old_saturation_solution, old_grad_saturation_solution_values); + saturation_fe_values.get_function_grads (old_old_saturation_solution, old_old_grad_saturation_solution_values); + darcy_fe_values.get_function_values (darcy_solution, present_darcy_solution_values); + + const double nu + = compute_viscosity (old_saturation_solution_values, + old_old_saturation_solution_values, + old_grad_saturation_solution_values, + old_old_grad_saturation_solution_values, + present_darcy_solution_values, + global_u_infty, + global_S_variation, + global_Omega_diameter, + saturation_fe_values.get_cell()->diameter(), + old_time_step, + viscosity); + + for (unsigned int q=0; q present_u; + for (unsigned int d=0; d grad_phi_i_s = saturation_fe_values.shape_grad (i, q); + + local_rhs(i) += (time_step * + f_saturation(old_s,viscosity) * + present_u * + grad_phi_i_s + - + time_step * + nu * + old_grad_saturation_solution_values[q] * grad_phi_i_s + + + old_s * phi_i_s) + * + saturation_fe_values.JxW(q); + } + + saturation_constraints.distribute_local_to_global (local_rhs, + local_dof_indices, + saturation_rhs); +} + +template +void +TwoPhaseFlowProblem:: +assemble_saturation_rhs_boundary_term (const FEFaceValues &saturation_fe_face_values, + const FEFaceValues &darcy_fe_face_values, + const std::vector &local_dof_indices) +{ + const unsigned int dofs_per_cell = saturation_fe_face_values.dofs_per_cell; + const unsigned int n_face_q_points = saturation_fe_face_values.n_quadrature_points; + + Vector local_rhs (dofs_per_cell); + + std::vector old_saturation_solution_values_face(n_face_q_points); + std::vector > present_darcy_solution_values_face(n_face_q_points, Vector(dim+1)); + std::vector neighbor_saturation (n_face_q_points); + + saturation_fe_face_values.get_function_values (old_saturation_solution, old_saturation_solution_values_face); + darcy_fe_face_values.get_function_values (darcy_solution, present_darcy_solution_values_face); + + SaturationBoundaryValues saturation_boundary_values; + saturation_boundary_values + .value_list (saturation_fe_face_values.get_quadrature_points(), + neighbor_saturation); + + for (unsigned int q=0; q present_u_face; + for (unsigned int d=0; d= 0); + + for (unsigned int i=0; i +void TwoPhaseFlowProblem::solve () +{ + solve_pressure_velocity_part = determine_whether_to_solve_pressure_velocity_part (); + + if ( timestep_number <= 3 || solve_pressure_velocity_part == true ) + { + std::cout << " Solving darcy system (pressure-velocity part)..." << std::endl; + + assemble_darcy_system (); + build_darcy_preconditioner (); + + { + const LinearSolvers::InverseMatrix + mp_inverse (darcy_preconditioner_matrix.block(1,1), *Mp_preconditioner); + + const LinearSolvers::BlockSchurPreconditioner + preconditioner (darcy_matrix, mp_inverse, *Amg_preconditioner); + + SolverControl solver_control (darcy_matrix.m(), + 1e-6*darcy_rhs.l2_norm()); + + SolverGMRES + gmres (solver_control, + SolverGMRES::AdditionalData(100)); + + for (unsigned int i=0; i cg (solver_control); + + TrilinosWrappers::PreconditionIC preconditioner; + preconditioner.initialize (saturation_matrix); + + cg.solve (saturation_matrix, saturation_solution, + saturation_rhs, preconditioner); + + + saturation_constraints.distribute (saturation_solution); + + project_back_saturation (); + + std::cout << " " + << solver_control.last_step() + << " CG iterations for saturation." + << std::endl; + + } + +} + +template +bool +TwoPhaseFlowProblem::determine_whether_to_solve_pressure_velocity_part () const +{ + if (timestep_number <= 3) + return true; + + const QGauss quadrature_formula(saturation_degree+2); + const unsigned int n_q_points = quadrature_formula.size(); + + FEValues fe_values (saturation_fe, quadrature_formula, + update_values | update_quadrature_points); + + std::vector old_saturation_after_solving_pressure (n_q_points); + std::vector present_saturation (n_q_points); + + const RandomMedium::KInverse k_inverse; +// const SingleCurvingCrack::KInverse k_inverse; + + std::vector > k_inverse_values (n_q_points); + + double max_global_aop_indicator = 0.0; + + typename DoFHandler::active_cell_iterator + cell = saturation_dof_handler.begin_active(), + endc = saturation_dof_handler.end(); + for (; cell!=endc; ++cell) + { + double max_local_mobility_reciprocal_difference = 0.0; + double max_local_permeability_inverse_l1_norm = 0.0; + + fe_values.reinit(cell); + fe_values.get_function_values (nth_saturation_solution_after_solving_pressure_part, + old_saturation_after_solving_pressure); + fe_values.get_function_values (saturation_solution, + present_saturation); + + k_inverse.value_list (fe_values.get_quadrature_points(), + k_inverse_values); + + for (unsigned int q=0; q 5.0 ) + { + return true; + } + else + { + std::cout << " Activating adaptive operating splitting" << std::endl; + return false; + } +} + +template +void +TwoPhaseFlowProblem:: +compute_refinement_indicators (Vector &refinement_indicators) const +{ + + const QMidpoint quadrature_formula; + FEValues fe_values (saturation_fe, quadrature_formula, update_gradients); + std::vector > grad_saturation (1); + + double max_refinement_indicator = 0.0; + + typename DoFHandler::active_cell_iterator + cell = saturation_dof_handler.begin_active(), + endc = saturation_dof_handler.end(); + for (unsigned int cell_no=0; cell!=endc; ++cell, ++cell_no) + { + fe_values.reinit(cell); + fe_values.get_function_grads (predictor_saturation_solution, + grad_saturation); + + refinement_indicators(cell_no) + = std::log( 1.0 + std::sqrt( grad_saturation[0] * + grad_saturation[0] ) ); + max_refinement_indicator = std::max(max_refinement_indicator, + refinement_indicators(cell_no)); + } + +// std::cout << "max_refinement_indicator =" << max_refinement_indicator << std::endl; +} + +template +void +TwoPhaseFlowProblem:: +refine_grid (const Vector &refinement_indicators) +{ + const double current_saturation_level = saturation_level + + n_refinement_steps; + + { + typename DoFHandler::active_cell_iterator + cell = saturation_dof_handler.begin_active(), + endc = saturation_dof_handler.end(); + + for (unsigned int cell_no=0; cell!=endc; ++cell, ++cell_no) + { + cell->clear_coarsen_flag(); + cell->clear_refine_flag(); + + if ((cell->level() < current_saturation_level) && + (std::fabs(refinement_indicators(cell_no)) > saturation_value)) + cell->set_refine_flag(); + else + if ((cell->level() > double(n_refinement_steps)) && + (std::fabs(refinement_indicators(cell_no)) < 0.75 * saturation_value)) + cell->set_coarsen_flag(); + } + } + + triangulation.prepare_coarsening_and_refinement (); + + unsigned int number_of_cells_refine = 0; + unsigned int number_of_cells_coarsen = 0; + + { + typename DoFHandler::active_cell_iterator + cell = saturation_dof_handler.begin_active(), + endc = saturation_dof_handler.end(); + + for (; cell!=endc; ++cell) + if (cell->refine_flag_set()) + ++number_of_cells_refine; + else + if (cell->coarsen_flag_set()) + ++number_of_cells_coarsen; + } + + std::cout << " " + << number_of_cells_refine + << " cell(s) are going to be refined." + << std::endl; + std::cout << " " + << number_of_cells_coarsen + << " cell(s) are going to be coarsened." + << std::endl; + + std::cout << std::endl; + + if ( number_of_cells_refine > 0 || number_of_cells_coarsen > 0 ) + { + std::vector x_saturation (3); + x_saturation[0] = saturation_solution; + x_saturation[1] = old_saturation_solution; + x_saturation[2] = nth_saturation_solution_after_solving_pressure_part; + + std::vector x_darcy (2); + x_darcy[0] = nth_darcy_solution_after_solving_pressure_part; + x_darcy[1] = n_minus_oneth_darcy_solution_after_solving_pressure_part; + + SolutionTransfer saturation_soltrans(saturation_dof_handler); + + SolutionTransfer darcy_soltrans(darcy_dof_handler); + + + triangulation.prepare_coarsening_and_refinement(); + saturation_soltrans.prepare_for_coarsening_and_refinement(x_saturation); + + darcy_soltrans.prepare_for_coarsening_and_refinement(x_darcy); + + triangulation.execute_coarsening_and_refinement (); + setup_dofs (); + + std::vector tmp_saturation (3); + tmp_saturation[0].reinit (saturation_solution); + tmp_saturation[1].reinit (saturation_solution); + tmp_saturation[2].reinit (saturation_solution); + saturation_soltrans.interpolate(x_saturation, tmp_saturation); + + saturation_solution = tmp_saturation[0]; + old_saturation_solution = tmp_saturation[1]; + nth_saturation_solution_after_solving_pressure_part = tmp_saturation[2]; + + std::vector tmp_darcy (2); + tmp_darcy[0].reinit (darcy_solution); + tmp_darcy[1].reinit (darcy_solution); + darcy_soltrans.interpolate(x_darcy, tmp_darcy); + + nth_darcy_solution_after_solving_pressure_part = tmp_darcy[0]; + n_minus_oneth_darcy_solution_after_solving_pressure_part = tmp_darcy[1]; + + rebuild_saturation_matrix = true; + } + else + { + rebuild_saturation_matrix = false; + + std::vector darcy_block_component (dim+1,0); + darcy_block_component[dim] = 1; + + std::vector darcy_dofs_per_block (2); + DoFTools::count_dofs_per_block (darcy_dof_handler, darcy_dofs_per_block, darcy_block_component); + const unsigned int n_u = darcy_dofs_per_block[0], + n_p = darcy_dofs_per_block[1], + n_s = saturation_dof_handler.n_dofs(); + + std::cout << "Number of active cells: " + << triangulation.n_active_cells() + << " (on " + << triangulation.n_levels() + << " levels)" + << std::endl + << "Number of degrees of freedom: " + << n_u + n_p + n_s + << " (" << n_u << '+' << n_p << '+'<< n_s <<')' + << std::endl + << std::endl; + } + +} + + +template +void TwoPhaseFlowProblem::output_results () const +{ + if ( solve_pressure_velocity_part == false ) + return; + + const FESystem joint_fe (darcy_fe, 1, + saturation_fe, 1); + DoFHandler joint_dof_handler (triangulation); + joint_dof_handler.distribute_dofs (joint_fe); + Assert (joint_dof_handler.n_dofs() == + darcy_dof_handler.n_dofs() + saturation_dof_handler.n_dofs(), + ExcInternalError()); + + Vector joint_solution (joint_dof_handler.n_dofs()); + + { + std::vector local_joint_dof_indices (joint_fe.dofs_per_cell); + std::vector local_darcy_dof_indices (darcy_fe.dofs_per_cell); + std::vector local_saturation_dof_indices (saturation_fe.dofs_per_cell); + + typename DoFHandler::active_cell_iterator + joint_cell = joint_dof_handler.begin_active(), + joint_endc = joint_dof_handler.end(), + darcy_cell = darcy_dof_handler.begin_active(), + saturation_cell = saturation_dof_handler.begin_active(); + + for (; joint_cell!=joint_endc; ++joint_cell, ++darcy_cell, ++saturation_cell) + { + joint_cell->get_dof_indices (local_joint_dof_indices); + darcy_cell->get_dof_indices (local_darcy_dof_indices); + saturation_cell->get_dof_indices (local_saturation_dof_indices); + + for (unsigned int i=0; i joint_solution_names; + switch (dim) + { + case 2: + joint_solution_names.push_back ("u"); + joint_solution_names.push_back ("v"); + break; + + case 3: + joint_solution_names.push_back ("u"); + joint_solution_names.push_back ("v"); + joint_solution_names.push_back ("w"); + break; + + default: + Assert (false, ExcNotImplemented()); + } + joint_solution_names.push_back ("pressure"); + joint_solution_names.push_back ("saturation"); + + std::vector + data_component_interpretation + (dim, DataComponentInterpretation::component_is_part_of_vector); + data_component_interpretation + .push_back (DataComponentInterpretation::component_is_scalar); + data_component_interpretation + .push_back (DataComponentInterpretation::component_is_scalar); + + DataOut data_out; + + data_out.attach_dof_handler (joint_dof_handler); + data_out.add_data_vector (joint_solution, joint_solution_names, + DataOut::type_dof_data, + data_component_interpretation); + + data_out.build_patches (); + + std::string filename = "solution-" + + Utilities::int_to_string (timestep_number, 5) + ".tec"; + std::ofstream output (filename.c_str()); + data_out.write_tecplot (output); +} + + +template +void +TwoPhaseFlowProblem::project_back_saturation () +{ + for (unsigned int i=0; i 1) + saturation_solution(i) = 1; +} + + +template +double +TwoPhaseFlowProblem::get_maximal_velocity () const +{ + QGauss quadrature_formula(darcy_degree+2); + const unsigned int n_q_points + = quadrature_formula.size(); + + FEValues darcy_fe_values (darcy_fe, quadrature_formula, + update_values); + std::vector > darcy_solution_values(n_q_points, + Vector(dim+1)); + double max_velocity = 0; + + typename DoFHandler::active_cell_iterator + cell = darcy_dof_handler.begin_active(), + endc = darcy_dof_handler.end(); + for (; cell!=endc; ++cell) + { + darcy_fe_values.reinit (cell); + darcy_fe_values.get_function_values (darcy_solution, darcy_solution_values); + + for (unsigned int q=0; q velocity; + for (unsigned int i=0; i +std::pair +TwoPhaseFlowProblem::get_extrapolated_saturation_range () const +{ + const QGauss quadrature_formula(saturation_degree+2); + const unsigned int n_q_points = quadrature_formula.size(); + + FEValues fe_values (saturation_fe, quadrature_formula, + update_values); + std::vector old_saturation_values(n_q_points); + std::vector old_old_saturation_values(n_q_points); + + if (timestep_number != 0) + { + double min_saturation = (1. + time_step/old_time_step) * + old_saturation_solution.linfty_norm() + + + time_step/old_time_step * + old_old_saturation_solution.linfty_norm(), + max_saturation = -min_saturation; + + typename DoFHandler::active_cell_iterator + cell = saturation_dof_handler.begin_active(), + endc = saturation_dof_handler.end(); + for (; cell!=endc; ++cell) + { + fe_values.reinit (cell); + fe_values.get_function_values (old_saturation_solution, + old_saturation_values); + fe_values.get_function_values (old_old_saturation_solution, + old_old_saturation_values); + + for (unsigned int q=0; q::active_cell_iterator + cell = saturation_dof_handler.begin_active(), + endc = saturation_dof_handler.end(); + for (; cell!=endc; ++cell) + { + fe_values.reinit (cell); + fe_values.get_function_values (old_saturation_solution, + old_saturation_values); + + for (unsigned int q=0; q +double +TwoPhaseFlowProblem:: +compute_viscosity (const std::vector &old_saturation, + const std::vector &old_old_saturation, + const std::vector > &old_saturation_grads, + const std::vector > &old_old_saturation_grads, + const std::vector > &present_darcy_values, + const double global_u_infty, + const double global_S_variation, + const double global_Omega_diameter, + const double cell_diameter, + const double old_time_step, + const double viscosity) +{ + const double beta = 0.08 * dim; + const double alpha = 1; + + if (global_u_infty == 0) + return 5e-3 * cell_diameter; + + const unsigned int n_q_points = old_saturation.size(); + + double max_residual = 0; + double max_velocity = 0; + + for (unsigned int q=0; q < n_q_points; ++q) + { + Tensor<1,dim> u; + for (unsigned int d=0; d +void TwoPhaseFlowProblem::run () +{ + unsigned int pre_refinement_step = 0; + + GridGenerator::hyper_cube (triangulation, 0, 1); + triangulation.refine_global (n_refinement_steps); + + setup_dofs (); + + start_time_iteration: + + VectorTools::project (saturation_dof_handler, + saturation_constraints, + QGauss(saturation_degree+2), + SaturationInitialValues(), + old_saturation_solution); + + timestep_number = 0; + double time = 0; + + do + { + std::cout << "Timestep " << timestep_number + << ": t=" << time + << ", dt=" << time_step + << std::endl; + + solve (); + + output_results (); + + solve_pressure_velocity_part = false; + + if ((timestep_number == 0) && + (pre_refinement_step < saturation_level)) + { + predictor_saturation_solution = saturation_solution; + predictor_saturation_solution.sadd (2.0, -1.0, old_saturation_solution); + Vector refinement_indicators (triangulation.n_active_cells()); + compute_refinement_indicators(refinement_indicators); + refine_grid(refinement_indicators); + ++pre_refinement_step; + goto start_time_iteration; + } + else + { + predictor_saturation_solution = saturation_solution; + predictor_saturation_solution.sadd (2.0, -1.0, old_saturation_solution); + Vector refinement_indicators (triangulation.n_active_cells()); + compute_refinement_indicators(refinement_indicators); + refine_grid(refinement_indicators); + } + + time += time_step; + ++timestep_number; + + old_old_saturation_solution = old_saturation_solution; + old_saturation_solution = saturation_solution; + + } + while (time <= 250); +} + + +int main () +{ + try + { + deallog.depth_console (0); + + TwoPhaseFlowProblem<3> two_phase_flow_problem(1); + two_phase_flow_problem.run (); + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + + return 0; +}