From: bangerth Date: Thu, 14 Aug 2008 13:11:40 +0000 (+0000) Subject: Add the variables we need for the BDF-2 time stepping scheme. Remove the code that... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=211b600b603e2da1387ce8b65d81c8a702d0000a;p=dealii-svn.git Add the variables we need for the BDF-2 time stepping scheme. Remove the code that dealt with inflow boundary conditions for the temperature field -- it didn't yield a contribution so far because the velocity was always tangential to the boundary, but we also can't enforce this sort of boundary condition any more because we now have diffusion and so we need Dirichlet or Neumann b.c. all around. git-svn-id: https://svn.dealii.org/trunk@16542 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-31/step-31.cc b/deal.II/examples/step-31/step-31.cc index 26e07b896c..3a452c2fd4 100644 --- a/deal.II/examples/step-31/step-31.cc +++ b/deal.II/examples/step-31/step-31.cc @@ -413,10 +413,12 @@ class BoussinesqFlowProblem BlockSparseMatrix preconditioner_matrix; double time_step; + double old_time_step; unsigned int timestep_number; BlockVector solution; BlockVector old_solution; + BlockVector old_old_solution; BlockVector system_rhs; //boost::shared_ptr::type> A_preconditioner; @@ -876,6 +878,8 @@ BoussinesqFlowProblem::BoussinesqFlowProblem (const unsigned int degree) FE_Q(degree), 1), dof_handler (triangulation), time_step (0), + old_time_step (0), + timestep_number (0), rebuild_matrices (true), rebuild_preconditioner (true) {} @@ -1207,6 +1211,12 @@ void BoussinesqFlowProblem::setup_dofs (const bool setup_matrices) old_solution.block(2).reinit (n_T); old_solution.collect_sizes (); + old_old_solution.reinit (3); + old_old_solution.block(0).reinit (n_u); + old_old_solution.block(1).reinit (n_p); + old_old_solution.block(2).reinit (n_T); + old_old_solution.collect_sizes (); + system_rhs.reinit (3); system_rhs.block(0).reinit (n_u); system_rhs.block(1).reinit (n_p); @@ -1762,23 +1772,16 @@ void BoussinesqFlowProblem::assemble_system () template void BoussinesqFlowProblem::assemble_rhs_T () { + // TODO: right now, always do explicit + // Euler + const bool is_first_timestep = (timestep_number == 0 ? true : true); + + QGauss quadrature_formula(degree+2); QGauss face_quadrature_formula(degree+2); FEValues fe_values (fe, quadrature_formula, update_values | update_gradients | update_quadrature_points | update_JxW_values); - FEFaceValues fe_face_values (fe, face_quadrature_formula, - update_values | update_normal_vectors | - update_quadrature_points | - update_JxW_values); - FESubfaceValues fe_subface_values (fe, face_quadrature_formula, - update_values | - update_normal_vectors | - update_JxW_values); - FEFaceValues fe_face_values_neighbor (fe, face_quadrature_formula, - update_values); - FESubfaceValues fe_subface_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.size(); @@ -1802,18 +1805,8 @@ void BoussinesqFlowProblem::assemble_rhs_T () // part of the FE system. 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 > > old_solution_grads( n_q_points, std::vector >(dim+2)); @@ -1870,45 +1863,6 @@ void BoussinesqFlowProblem::assemble_rhs_T () fe_values.JxW(q); } - for (unsigned int face_no=0; face_no::faces_per_cell; - ++face_no) - if (cell->at_boundary(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); - - temperature_boundary_values - .value_list (fe_face_values.get_quadrature_points(), - neighbor_temperature); - - for (unsigned int q=0; q present_u_face; - for (unsigned int d=0; d= 0); - - for (unsigned int i=0; iget_dof_indices (local_dof_indices); for (unsigned int i=0; i::solve () solution.block(0) = up.block(0); solution.block(1) = up.block(1); } - // for DGQ1 needs to be /15 + + // TODO: determine limit of stability for + // the time step (whether it needs to be /4 + // or whether we could get away with a + // bigger time step) + old_time_step = time_step; time_step = GridTools::minimal_cell_diameter(triangulation) / std::max (get_maximal_velocity(), .05) / 4; @@ -2000,15 +1959,8 @@ void BoussinesqFlowProblem::solve () PreconditionSSOR<> preconditioner; preconditioner.initialize (system_matrix.block(2,2), 1.2); - try - { - cg.solve (system_matrix.block(2,2), solution.block(2), - system_rhs.block(2), preconditioner); - } - catch (...) - { - abort (); - } + cg.solve (system_matrix.block(2,2), solution.block(2), + system_rhs.block(2), preconditioner); // produce a consistent temperature field hanging_node_constraints.distribute (solution); @@ -2094,6 +2046,8 @@ void BoussinesqFlowProblem::refine_mesh () Vector x_old_solution (dof_handler.n_dofs()); x_old_solution = old_solution; + Assert (false, ExcMessage ("Need to do the same for old_old_solution")); + soltrans.prepare_for_coarsening_and_refinement(x_old_solution); triangulation.execute_coarsening_and_refinement (); @@ -2241,7 +2195,8 @@ void BoussinesqFlowProblem::run () time += time_step; ++timestep_number; - old_solution = solution; + old_old_solution = old_solution; + old_solution = solution; std::cout << std::endl;