From f9775b6b2beb452a5b75af120d271c6516c68396 Mon Sep 17 00:00:00 2001 From: bangerth Date: Wed, 10 Oct 2007 02:17:43 +0000 Subject: [PATCH] Get the thing to run. Simplify a bunch of places by removing stuff that was specific to step-21. git-svn-id: https://svn.dealii.org/trunk@15276 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-22/step-22.cc | 241 +++++----------------------- 1 file changed, 41 insertions(+), 200 deletions(-) diff --git a/deal.II/examples/step-22/step-22.cc b/deal.II/examples/step-22/step-22.cc index d6cc6b48db..ae8f6d8275 100644 --- a/deal.II/examples/step-22/step-22.cc +++ b/deal.II/examples/step-22/step-22.cc @@ -83,7 +83,6 @@ class BoussinesqFlowProblem double time_step; unsigned int timestep_number; - double viscosity; BlockVector solution; BlockVector old_solution; @@ -220,139 +219,6 @@ InitialValues::vector_value (const Point &p, - -namespace SingleCurvingCrack -{ - template - class KInverse : public TensorFunction<2,dim> - { - public: - 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 Tensor<1,dim> extract_u (const FEValuesBase &fe_values, @@ -470,8 +336,16 @@ void InverseMatrix::vmult (Vector &dst, SolverCG<> cg (solver_control, vector_memory); dst = 0; - - cg.solve (*matrix, dst, src, PreconditionIdentity()); + + try + { + cg.solve (*matrix, dst, src, PreconditionIdentity()); + } + catch (...) + { + std::cout << "ups 1..." << std::endl; + Assert (false, ExcInternalError()); + } } @@ -515,44 +389,6 @@ void SchurComplement::vmult (Vector &dst, -class ApproximateSchurComplement : public Subscriptor -{ - public: - ApproximateSchurComplement (const BlockSparseMatrix &A); - - void vmult (Vector &dst, - const Vector &src) const; - - private: - const SmartPointer > system_matrix; - - mutable Vector tmp1, tmp2; -}; - - -ApproximateSchurComplement:: -ApproximateSchurComplement (const BlockSparseMatrix &A) - : - system_matrix (&A), - tmp1 (A.block(0,0).m()), - tmp2 (A.block(0,0).m()) -{} - - -void ApproximateSchurComplement::vmult (Vector &dst, - const Vector &src) const -{ - system_matrix->block(0,1).vmult (tmp1, src); - system_matrix->block(0,0).precondition_Jacobi (tmp2, tmp1); - system_matrix->block(1,0).vmult (dst, tmp2); -} - - - - - - - template BoussinesqFlowProblem::BoussinesqFlowProblem (const unsigned int degree) : @@ -561,9 +397,8 @@ BoussinesqFlowProblem::BoussinesqFlowProblem (const unsigned int degree) FE_DGQ(degree), 1, FE_DGQ(degree), 1), dof_handler (triangulation), - n_refinement_steps (5), - time_step (0), - viscosity (0.2) + n_refinement_steps (3), + time_step (0) {} @@ -679,13 +514,11 @@ void BoussinesqFlowProblem::assemble_system () const PressureRightHandSide pressure_right_hand_side; const Buoyancy buoyancy; const PressureBoundaryValues pressure_boundary_values; - const RandomMedium::KInverse k_inverse; std::vector pressure_rhs_values (n_q_points); std::vector > buoyancy_values (n_q_points, Vector(dim)); std::vector boundary_values (n_face_q_points); - std::vector > k_inverse_values (n_q_points); std::vector > old_solution_values(n_q_points, Vector(dim+2)); std::vector > > old_solution_grads(n_q_points, @@ -706,9 +539,6 @@ void BoussinesqFlowProblem::assemble_system () pressure_rhs_values); buoyancy.vector_value_list (fe_values.get_quadrature_points(), buoyancy_values); - k_inverse.value_list (fe_values.get_quadrature_points(), - k_inverse_values); - for (unsigned int q=0; q::assemble_rhs_S () const Tensor<1,dim> grad_phi_i_s = extract_grad_s(fe_values, i, q); local_rhs(i) += (time_step * - f_saturation(old_s,viscosity) * + old_s * present_u * grad_phi_i_s + @@ -905,12 +735,11 @@ void BoussinesqFlowProblem::assemble_rhs_S () for (unsigned int i=0; i::solve () SchurComplement schur_complement (system_matrix, m_inverse); - ApproximateSchurComplement - approximate_schur_complement (system_matrix); - - InverseMatrix - preconditioner (approximate_schur_complement); - - SolverControl solver_control (system_matrix.block(0,0).m(), 1e-12*schur_rhs.l2_norm()); SolverCG<> cg (solver_control); - cg.solve (schur_complement, solution.block(1), schur_rhs, - preconditioner); + try + { + cg.solve (schur_complement, solution.block(1), schur_rhs, + PreconditionIdentity()); + } + catch (...) + { + std::cout << "ups 2..." << std::endl; + abort (); + } + std::cout << " " << solver_control.last_step() @@ -981,8 +812,18 @@ void BoussinesqFlowProblem::solve () SolverControl solver_control (system_matrix.block(2,2).m(), 1e-8*system_rhs.block(2).l2_norm()); SolverCG<> cg (solver_control); - cg.solve (system_matrix.block(2,2), solution.block(2), system_rhs.block(2), - PreconditionIdentity()); + + try + { + cg.solve (system_matrix.block(2,2), solution.block(2), system_rhs.block(2), + PreconditionIdentity()); + } + catch (...) + { + std::cout << "ups 3..." << std::endl; + abort (); + } + project_back_saturation (); -- 2.39.5