From c07dd330d8224d904718de27b1679e93a602ebf7 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Tue, 18 Aug 2009 12:06:37 +0000 Subject: [PATCH] Use Dirichlet boundary conditions for temperature instead of natural conditions. Use ConstraintMatrix for implementing this. Comment the things that have changed. git-svn-id: https://svn.dealii.org/trunk@19300 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-32/step-32.cc | 186 +++++++++++++++++++++------- 1 file changed, 144 insertions(+), 42 deletions(-) diff --git a/deal.II/examples/step-32/step-32.cc b/deal.II/examples/step-32/step-32.cc index b54dbe5a9b..e757edb002 100644 --- a/deal.II/examples/step-32/step-32.cc +++ b/deal.II/examples/step-32/step-32.cc @@ -138,8 +138,8 @@ namespace EquationData const double rho = (r-R0)/h; - Assert (rho >= -.05, ExcInternalError()); - Assert (rho <= 1, ExcInternalError()); + Assert (rho >= -.001, ExcInternalError()); + Assert (rho <= 1.001, ExcInternalError()); return T1+(T0-T1)*((1-rho)*(1-rho)); } @@ -680,6 +680,7 @@ namespace Assembly Vector local_rhs; std::vector local_dof_indices; + FullMatrix matrix_for_bc; }; template @@ -687,7 +688,9 @@ namespace Assembly TemperatureRHS (const FiniteElement &temperature_fe) : local_rhs (temperature_fe.dofs_per_cell), - local_dof_indices (temperature_fe.dofs_per_cell) + local_dof_indices (temperature_fe.dofs_per_cell), + matrix_for_bc (temperature_fe.dofs_per_cell, + temperature_fe.dofs_per_cell) {} @@ -696,7 +699,8 @@ namespace Assembly TemperatureRHS (const TemperatureRHS &data) : local_rhs (data.local_rhs), - local_dof_indices (data.local_dof_indices) + local_dof_indices (data.local_dof_indices), + matrix_for_bc (data.matrix_for_bc) {} } } @@ -1302,6 +1306,47 @@ compute_viscosity (const std::vector &old_temperature, // side is cheap anyway so we won't even // notice that this part is not parallized // by threads. + // + // Regarding the implementation of + // inhomogeneous Dirichlet boundary + // conditions: Since we use the temperature + // ConstraintMatrix, we can apply the + // boundary conditions directly when + // building the respective matrix and right + // hand side. In this case, the boundary + // conditions are inhomogeneous, which + // makes this procedure somewhat + // tricky. Remember that we get the matrix + // from some other function. However, the + // correct imposition of boundary + // conditions needs the matrix data we work + // on plus the right hand side + // simultaneously, since the right hand + // side is created by Gaussian elimination + // on the matrix rows. In order to not + // introduce the matrix assembly at this + // place, but still having the matrix data + // available, we choose to create a dummy + // matrix matrix_for_bc that + // we only fill with data when we need it + // for imposing boundary conditions. These + // positions are exactly those where we + // have an inhomogeneous entry in the + // ConstraintMatrix. There are only a few + // such positions (on the boundary dofs), + // so it is still much cheaper to use this + // function than to create the full matrix + // here. To implement this, we ask the + // constraint matrix whether the dof under + // consideration is inhomogeneously + // constraint. In that case, we generate + // the respective matrix column that we + // need for creating the correct right hand + // side. Note that this (manually + // generated) matrix entry needs to be + // exactly the entry that we would fill the + // matrix with — otherwise, this will + // not work. template void BoussinesqFlowProblem::project_temperature_field () { @@ -1318,6 +1363,7 @@ void BoussinesqFlowProblem::project_temperature_field () std::vector local_dof_indices (dofs_per_cell); Vector cell_vector (dofs_per_cell); + FullMatrix matrix_for_bc (dofs_per_cell, dofs_per_cell); typename DoFHandler::active_cell_iterator cell = temperature_dof_handler.begin_active(), @@ -1333,34 +1379,42 @@ void BoussinesqFlowProblem::project_temperature_field () if (cell->subdomain_id() == Utilities::Trilinos::get_this_mpi_process(trilinos_communicator)) { + cell->get_dof_indices (local_dof_indices); fe_values.reinit(cell); - const std::vector &weights = fe_values.get_JxW_values (); EquationData::TemperatureInitialValues().value_list (fe_values.get_quadrature_points(), rhs_values); cell_vector = 0; + matrix_for_bc = 0; for (unsigned int point=0; pointget_dof_indices (local_dof_indices); + { + cell_vector(i) += rhs_values[point] * + fe_values.shape_value(i,point) * + fe_values.JxW(point); + if (temperature_constraints.is_inhomogeneously_constrained(local_dof_indices[i])) + { + for (unsigned int j=0; j memory; - SolverCG cg(control,memory); + SolverControl solver_control(5*rhs.size(), 1e-12*rhs.l2_norm()); + SolverCG cg(solver_control); - TrilinosWrappers::PreconditionIC preconditioner_mass; - preconditioner_mass.initialize(temperature_mass_matrix); + TrilinosWrappers::PreconditionSSOR preconditioner_mass; + preconditioner_mass.initialize(temperature_mass_matrix, 1.3); cg.solve (temperature_mass_matrix, solution, rhs, preconditioner_mass); @@ -1559,7 +1613,9 @@ void BoussinesqFlowProblem::setup_temperature_matrices () // that part of the matrix that is roughly // equal to the degrees of freedom located // on those cells that it will actually - // work on. + // work on. Note how we set boundary + // conditions on the temperature by using + // the ConstraintMatrix object. // // After this, we have to set up the // various partitioners (of type @@ -1641,6 +1697,14 @@ void BoussinesqFlowProblem::setup_dofs () DoFRenumbering::subdomain_wise (temperature_dof_handler); temperature_constraints.clear (); + VectorTools::interpolate_boundary_values (temperature_dof_handler, + 0, + EquationData::TemperatureInitialValues(), + temperature_constraints); + VectorTools::interpolate_boundary_values (temperature_dof_handler, + 1, + EquationData::TemperatureInitialValues(), + temperature_constraints); DoFTools::make_hanging_node_constraints (temperature_dof_handler, temperature_constraints); temperature_constraints.close (); @@ -2271,7 +2335,19 @@ void BoussinesqFlowProblem::assemble_temperature_matrix () // points (which are necessary for // calculating the artificial viscosity of // stabilization), but is otherwise similar - // to the other assembly functions. + // to the other assembly functions. Notice, + // once again, how we resolve the dilemma + // of having inhomogeneous boundary + // conditions, but just making a right hand + // side at this point (compare the comments + // for the project function): We create + // some matrix columns with exactly the + // values that would be entered for the + // temperature stiffness matrix, in case we + // have inhomogeneously constrained + // dofs. That will account for the correct + // balance of the right hand side vector + // with the matrix system of temperature. template void BoussinesqFlowProblem:: local_assemble_temperature_rhs (const std::pair global_T_range, @@ -2290,6 +2366,8 @@ local_assemble_temperature_rhs (const std::pair global_T_range, const FEValuesExtractors::Vector velocities (0); data.local_rhs = 0; + data.matrix_for_bc = 0; + cell->get_dof_indices (data.local_dof_indices); scratch.temperature_fe_values.reinit (cell); @@ -2375,21 +2453,37 @@ local_assemble_temperature_rhs (const std::pair global_T_range, scratch.old_velocity_values[q]); for (unsigned int i=0; iget_dof_indices (data.local_dof_indices); } @@ -2400,7 +2494,8 @@ copy_local_to_global_temperature_rhs (const Assembly::CopyData::TemperatureRHS::solve () if (stokes_constraints.is_constrained (i)) distributed_stokes_solution(i) = 0; - SolverControl solver_control (stokes_matrix.m(), 1e-21*stokes_rhs.l2_norm()); SolverBicgstab bicgstab (solver_control, false); @@ -2592,10 +2686,16 @@ void BoussinesqFlowProblem::solve () old_time_step = time_step; const double maximal_velocity = get_maximal_velocity(); - time_step = 1./(1.8*dim*std::sqrt(1.*dim)) / - temperature_degree * - GridTools::minimal_cell_diameter(triangulation) / - maximal_velocity; + if (maximal_velocity > 1e-10) + time_step = 1./(1.8*dim*std::sqrt(1.*dim)) / + temperature_degree * + GridTools::minimal_cell_diameter(triangulation) / + maximal_velocity; + else + time_step = 1./(1.8*dim*std::sqrt(1.*dim)) / + temperature_degree * + GridTools::minimal_cell_diameter(triangulation) / + 1e-10; pcout << " " << "Time step: " << time_step/EquationData::year_in_seconds @@ -2612,7 +2712,7 @@ void BoussinesqFlowProblem::solve () { SolverControl solver_control (temperature_matrix.m(), - 1e-8*temperature_rhs.l2_norm()); + 1e-12*temperature_rhs.l2_norm()); SolverCG cg (solver_control); TrilinosWrappers::MPI::Vector @@ -2811,7 +2911,7 @@ void BoussinesqFlowProblem::refine_mesh (const unsigned int max_grid_level) GridRefinement::refine_and_coarsen_fixed_fraction (triangulation, estimated_error_per_cell, - 0.8, 0.1); + 0.6, 0.2); if (triangulation.n_levels() > max_grid_level) for (typename Triangulation::active_cell_iterator cell = triangulation.begin_active(max_grid_level); @@ -2846,6 +2946,8 @@ void BoussinesqFlowProblem::refine_mesh (const unsigned int max_grid_level) temperature_solution = tmp[0]; old_temperature_solution = tmp[1]; + temperature_constraints.distribute(temperature_solution); + temperature_constraints.distribute(old_temperature_solution); TrilinosWrappers::BlockVector x_stokes_new = stokes_solution; stokes_trans.interpolate (x_stokes, x_stokes_new); @@ -2919,8 +3021,6 @@ void BoussinesqFlowProblem::run () solve (); - output_results (); - pcout << std::endl; if ((timestep_number == 0) && @@ -2934,6 +3034,8 @@ void BoussinesqFlowProblem::run () if ((timestep_number > 0) && (timestep_number % 5 == 0)) refine_mesh (initial_refinement + n_pre_refinement_steps); + output_results (); + time += time_step; ++timestep_number; -- 2.39.5