From: Abner J. Salgado Date: Fri, 18 Sep 2009 23:35:57 +0000 (+0000) Subject: Turns out I was using wrong boundary conditions + wrong weak X-Git-Tag: v8.0.0~7097 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=f9aaf1d1675ea4aa1c423c9aab07f27aeae8b142;p=dealii.git Turns out I was using wrong boundary conditions + wrong weak formulation. Everything is fixed now and if we run the program it gives the correct results. Yay!!! git-svn-id: https://svn.dealii.org/trunk@19478 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-35/step-35.cc b/deal.II/examples/step-35/step-35.cc index ec48e880ac..ebbe40b196 100644 --- a/deal.II/examples/step-35/step-35.cc +++ b/deal.II/examples/step-35/step-35.cc @@ -347,10 +347,10 @@ namespace EquationData { const double Um = 1.5; const double H = 4.1; - return 4.*Um*p(1)* (H - p(1))/(H*H); + return 4.*Um*p(1)*(H - p(1))/(H*H); } else - return 0; + return 0.; } @@ -380,7 +380,7 @@ namespace EquationData double Pressure::value (const Point &p, const unsigned int) const { - return 0.; + return 25.-p(0); } template @@ -434,16 +434,15 @@ class NavierStokesProjection SparsityPattern sparsity_pattern_velocity, sparsity_pattern_pressure, sparsity_pattern_pres_vel; SparseMatrix vel_Laplace_plus_Mass, vel_it_matrix[dim], vel_Mass, vel_Laplace, vel_Advection, - pres_Laplace, pres_Mass, pres_Diff[dim]; - ConstraintMatrix pres_regularization; + pres_Laplace, pres_Mass, pres_Diff[dim], pres_iterative; Vector pres_n, pres_n_minus_1, phi_n, phi_n_minus_1, u_n[dim], u_n_minus_1[dim], u_star[dim], force[dim], v_tmp, pres_tmp, rot_u; - SparseILU prec_velocity[dim]; - SparseDirectUMFPACK prec_mass, prec_pressure, prec_vel_mass; + SparseILU prec_velocity[dim], prec_pres_Laplace; + SparseDirectUMFPACK prec_mass, prec_vel_mass; DeclException2 (ExcInvalidTimeStep, double, double, << " The time step " << arg1 << " is out of range." << std::endl << " The permitted range is (0," << arg2 << "]"); @@ -661,9 +660,9 @@ void NavierStokesProjection::Initialize() { vel_exact.set_time (t_0); vel_exact.set_component(d); - VectorTools::interpolate (dof_handler_velocity, vel_exact, u_n_minus_1[d]); + VectorTools::interpolate (dof_handler_velocity, ZeroFunction(), u_n_minus_1[d]); vel_exact.advance_time (dt); - VectorTools::interpolate (dof_handler_velocity, vel_exact, u_n[d]); + VectorTools::interpolate (dof_handler_velocity, ZeroFunction(), u_n[d]); } } @@ -674,7 +673,7 @@ void NavierStokesProjection::Initialize() // initialize the sparsity patterns, // the constraints (if any) and // assemble the matrices that do not - // depend on the timestep $\Delta t$. + // depend on the timestep dt. template void NavierStokesProjection::init_velocity_matrices() { @@ -716,20 +715,14 @@ void NavierStokesProjection::init_pressure_matrices() dof_handler_pressure.max_couplings_between_dofs()); DoFTools::make_sparsity_pattern (dof_handler_pressure, sparsity_pattern_pressure); - pres_regularization.clear(); - pres_regularization.add_line(0); - pres_regularization.close(); - pres_regularization.condense (sparsity_pattern_pressure); - sparsity_pattern_pressure.compress(); pres_Laplace.reinit (sparsity_pattern_pressure); + pres_iterative.reinit (sparsity_pattern_pressure); pres_Mass.reinit (sparsity_pattern_pressure); MatrixCreator::create_laplace_matrix (dof_handler_pressure, quadrature_pressure, pres_Laplace); MatrixCreator::create_mass_matrix ( dof_handler_pressure, quadrature_pressure, pres_Mass); - - pres_regularization.condense (pres_Laplace); } @@ -757,7 +750,7 @@ void NavierStokesProjection::init_gradient_operator() InitGradPerTaskData per_task_data (0, fe_velocity.dofs_per_cell, fe_pressure.dofs_per_cell); InitGradScratchData scratch_data (fe_velocity, fe_pressure, quadrature_velocity, - update_values | update_JxW_values, update_gradients); + update_gradients | update_JxW_values, update_values); for (unsigned int d=0; d::assemble_one_cell_of_gradient (const SIterator { for (unsigned int i=0; i::copy_advection_local_to_global( template void NavierStokesProjection::projection_step (const bool reinit_prec) { - if (reinit_prec) - prec_pressure.initialize (pres_Laplace); - + pres_iterative.copy_from (pres_Laplace); + pres_tmp = 0.; for (unsigned d=0; d bval; + if (reinit_prec) + VectorTools::interpolate_boundary_values (dof_handler_pressure, 3, + ZeroFunction(), bval); + + MatrixTools::apply_boundary_values (bval, pres_iterative, phi_n, pres_tmp); + + if (reinit_prec) + prec_pres_Laplace.initialize(pres_iterative, + SparseILU::AdditionalData (vel_diag_strength, vel_off_diagonals) ); + + SolverControl solvercontrol (vel_max_its, vel_eps*pres_tmp.l2_norm()); + SolverCG<> cg (solvercontrol); + cg.solve (pres_iterative, phi_n, pres_tmp, prec_pres_Laplace); + + phi_n *= 1.5/dt; } @@ -1074,7 +1077,7 @@ void NavierStokesProjection::update_pressure (const bool reinit_prec) // step-31 and so I will not // elaborate on it. There is one // small detail here. It is often - // interested to see the vorticity of + // interesting to see the vorticity of // the flow. But, since we are using // it here only for plotting // purposes, we are not going to