From f1ece58ecbbc641d14d04896a1d9c1d1afc55130 Mon Sep 17 00:00:00 2001 From: "Abner J. Salgado" Date: Tue, 8 Sep 2009 17:11:54 +0000 Subject: [PATCH] Everything is computing fine. I will add the vorticity to the plot method and we should be good to go... git-svn-id: https://svn.dealii.org/trunk@19418 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-35/Makefile | 2 +- deal.II/examples/step-35/doc/intro.dox | 7 +- deal.II/examples/step-35/parameter-file.prm | 45 +-- deal.II/examples/step-35/step-35.cc | 351 +++++--------------- 4 files changed, 105 insertions(+), 300 deletions(-) diff --git a/deal.II/examples/step-35/Makefile b/deal.II/examples/step-35/Makefile index 5df8a9405b..d148efb30b 100644 --- a/deal.II/examples/step-35/Makefile +++ b/deal.II/examples/step-35/Makefile @@ -14,7 +14,7 @@ target = step-35 # run-time checking of parameters and internal states is performed, so # you should set this value to `on' while you develop your program, # and to `off' when running production computations. -debug-mode = on +debug-mode = off # As third field, we need to give the path to the top-level deal.II diff --git a/deal.II/examples/step-35/doc/intro.dox b/deal.II/examples/step-35/doc/intro.dox index 68930b4586..23cbe2b0ea 100644 --- a/deal.II/examples/step-35/doc/intro.dox +++ b/deal.II/examples/step-35/doc/intro.dox @@ -1,3 +1,8 @@ + This program grew out of a student project by Abner Salgado at Texas A&M University. Most of the work for this program is by him. + + +

Introduction

-

Introduction

+The purpose of this program is to show how to effectively solve the incompressible time-dependent Navier-Stokes equations. +In previous tutorial programs (see for instance ****) we have seen how to solve the time-independent Stokes equations using a Schur-complement approach. However \ No newline at end of file diff --git a/deal.II/examples/step-35/parameter-file.prm b/deal.II/examples/step-35/parameter-file.prm index 52cc5cee06..12eb9d3cde 100644 --- a/deal.II/examples/step-35/parameter-file.prm +++ b/deal.II/examples/step-35/parameter-file.prm @@ -5,18 +5,15 @@ set Method_Form = rotational subsection Physical data # In this subsection we declare the physical data # The initial and final time, and the Reynolds number - set initial_time = 0.0 - set final_time = 1.0 - set Reynolds = 1. + set initial_time = 0. + set final_time = 8. + set Reynolds = 100 end subsection Time step data - # In this subsection we declare the data that is to be used for the convergence tests - # To test convergence in time we run the algorithm with the time step set from - # 'initial_dt' to 'final_dt', each time dividing the current time step by 'dt_decrement' - set initial_dt = .1 - set final_dt = 5e-3 - set dt_decrement = 2. + # In this subsection we declare the data that is to be used for time discretization, + # i.e. the time step dt + set dt = 5e-3 end subsection Space discretization @@ -31,31 +28,15 @@ subsection Data solve velocity # In this section we declare the parameters that are going to control the solution process # for the velocity. set max_iterations = 1000 # maximal number of iterations that GMRES must make - set eps = 1e-12 # stopping criterion - set Krylov_size = 10 # size of the Krylov subspace to be used in GMRES + set eps = 1e-6 # stopping criterion + set Krylov_size = 30 # size of the Krylov subspace to be used in GMRES set off_diagonals = 60 # number of off diagonals that ILU must compute - set diag_strength = 0 # diagonal strengthening value - set update_prec = 20 # this number indicates how often the preconditioner must be updated -end - -subsection Data solve projection - # In this section we declare the parameters that are going to control the solution process - # for the projection step. - set max_iterations = 1000 # maximal number of iterations that CG must make - set eps = 1e-12 # stopping criterion - set off_diagonals = 100 # number of off diagonals that ILU must compute. This one can be big as it is computed only once - set diag_strength = 1.5 # diagonal strengthening value -end - -subsection Data solve pressure update - # In this section we declare the parameters that are going to control the solution process - # for the pressure update step in the case of a rotational form of the method - set max_iterations = 1000 # maximal number of iterations that CG must make - set eps = 1e-12 # stopping criterion - set off_diagonals = 100 # number of off diagonals that ILU must compute. This one can be big as it is computed only once - set diag_strength = 1e-5 # diagonal strengthening value + set diag_strength = 0.01 # diagonal strengthening value + set update_prec = 10 # this number indicates how often the preconditioner must be updated end +#The output frequency +set output = 1 #Finally we set the verbosity level -set verbose = false \ No newline at end of file +set verbose = true \ No newline at end of file diff --git a/deal.II/examples/step-35/step-35.cc b/deal.II/examples/step-35/step-35.cc index 6cf60dadb9..d440d910c9 100644 --- a/deal.II/examples/step-35/step-35.cc +++ b/deal.II/examples/step-35/step-35.cc @@ -21,7 +21,6 @@ #include #include #include -#include #include #include #include @@ -87,30 +86,19 @@ namespace RunTimeParameters{ void read_data( const char *filename ); Method_Formulation form; double initial_time, - final_time, - Reynolds; - double initial_dt, - final_dt, - dt_decrement; + final_time, + Reynolds; + double dt; unsigned int n_of_global_refines, - pressure_degree; + pressure_degree; unsigned int vel_max_iterations, - vel_Krylov_size, - vel_off_diagonals, - vel_update_prec; + vel_Krylov_size, + vel_off_diagonals, + vel_update_prec; double vel_eps, - vel_diag_strength; - unsigned int proj_max_iterations, - proj_off_diagonals; - double proj_eps, - proj_diag_strength; - unsigned int pres_max_iterations, - pres_off_diagonals; - double pres_eps, - pres_diag_strength; + vel_diag_strength; bool verbose; unsigned int output; - protected: ParameterHandler prm; }; @@ -118,7 +106,6 @@ namespace RunTimeParameters{ // In the constructor of this class we declare all the parameters. // The details of how this works have been discussed somewhere else *** // so let's not elaborate on that - Data_Storage::Data_Storage(){ prm.declare_entry( "Method_Form", "rotational", Patterns::Selection( "rotational|standard" ), " Used to select the type of method that we are going to use. " ); @@ -129,14 +116,11 @@ namespace RunTimeParameters{ prm.leave_subsection(); prm.enter_subsection( "Time step data" ); - prm.declare_entry( "initial_dt", "0.1", Patterns::Double( 0. ), " The initial time step size. " ); - prm.declare_entry( "final_dt", "5e-4", Patterns::Double( 0. ), " The final time step size. " ); - prm.declare_entry( "dt_decrement", "2.", Patterns::Double( 1.5 ), - " The factor by which the time step will be divided. " ); + prm.declare_entry( "dt", "5e-4", Patterns::Double( 0. ), " The time step size. " ); prm.leave_subsection(); prm.enter_subsection( "Space discretization" ); - prm.declare_entry( "n_of_refines", "5", Patterns::Integer( 1, 15), + prm.declare_entry( "n_of_refines", "0", Patterns::Integer( 0, 15), " The number of global refines we do on the mesh. " ); prm.declare_entry( "pressure_fe_degree", "1", Patterns::Integer( 1, 5 ), " The polynomial degree for the pressure space. " ); @@ -155,28 +139,10 @@ namespace RunTimeParameters{ " This number indicates how often we need to update the preconditioner" ); prm.leave_subsection(); - prm.enter_subsection( "Data solve projection" ); - prm.declare_entry( "max_iterations", "1000", Patterns::Integer( 1, 1000 ), - " The maximal number of iterations CG must make. " ); - prm.declare_entry( "eps", "1e-12", Patterns::Double( 0. ), " The stopping criterion. " ); - prm.declare_entry( "off_diagonals", "100", Patterns::Integer(1), - " The number of off-diagonal elements ILU must compute" ); - prm.declare_entry( "diag_strength", "0.1", Patterns::Double( 0. ), " Diagonal strengthening coefficient. " ); - prm.leave_subsection(); - - prm.enter_subsection( "Data solve pressure update" ); - prm.declare_entry( "max_iterations", "1000", Patterns::Integer( 1, 1000 ), - " The maximal number of iterations CG must make. " ); - prm.declare_entry( "eps", "1e-12", Patterns::Double( 0. ), " The stopping criterion. " ); - prm.declare_entry( "off_diagonals", "10", Patterns::Integer(0), - " The number of off-diagonal elements that ILU must compute" ); - prm.declare_entry( "diag_strength", "0.", Patterns::Double(0), " Diagonal strengthening coefficient" ); - prm.leave_subsection(); - prm.declare_entry( "verbose", "true", Patterns::Bool(), " This indicates whether the output of the solution process should be verbose. " ); - prm.declare_entry( "output", "10", Patterns::Integer(1), + prm.declare_entry( "output", "1", Patterns::Integer(1), " This indicates between how many time steps we print the solution. " ); } @@ -201,9 +167,7 @@ namespace RunTimeParameters{ prm.leave_subsection(); prm.enter_subsection( "Time step data" ); - initial_dt = prm.get_double( "initial_dt" ); - final_dt = prm.get_double( "final_dt" ); - dt_decrement = prm.get_double( "dt_decrement" ); + dt = prm.get_double( "dt" ); prm.leave_subsection(); prm.enter_subsection( "Space discretization" ); @@ -220,20 +184,6 @@ namespace RunTimeParameters{ vel_update_prec = prm.get_integer( "update_prec" ); prm.leave_subsection(); - prm.enter_subsection( "Data solve projection" ); - proj_max_iterations = prm.get_integer( "max_iterations" ); - proj_eps = prm.get_double( "eps" ); - proj_off_diagonals = prm.get_integer( "off_diagonals" ); - proj_diag_strength = prm.get_double( "diag_strength" ); - prm.leave_subsection(); - - prm.enter_subsection( "Data solve pressure update" ); - pres_max_iterations = prm.get_integer( "max_iterations" ); - pres_eps = prm.get_double( "eps" ); - pres_off_diagonals = prm.get_integer( "off_diagonals" ); - pres_diag_strength = prm.get_double( "diag_strength" ); - prm.leave_subsection(); - verbose = prm.get_bool( "verbose" ); output = prm.get_integer( "output" ); @@ -273,12 +223,8 @@ namespace EquationData{ public: Velocity( const double initial_time = 0.0 ); virtual double value( const Point &p, const unsigned int component = 0 ) const; - virtual Tensor<1,dim> gradient( const Point &p, const unsigned int component = 0 ) const; virtual void value_list( const std::vector< Point > &points, std::vector &values, const unsigned int component = 0 ) const; - virtual void gradient_list( const std::vector< Point > &points, - std::vector< Tensor<1,dim> > &gradients, - const unsigned int component = 0 ) const; }; template Velocity::Velocity( const double initial_time ): @@ -293,66 +239,25 @@ namespace EquationData{ } template inline double Velocity::value( const Point &p, const unsigned int ) const{ - double return_value = std::cos( Function::get_time() ); - switch( MultiComponentFunction::comp ){ - case 0: - return_value *= -p(1); - break; - case 1: - return_value *= p(0); - break; - default: - Assert( false, ExcNotImplemented() ); - }; + double return_value = 0., dist = std::sqrt( p.square() ); + static const double Um = 1.5, H = 4.1; + if( MultiComponentFunction::comp == 0 ) + return_value = 4.*Um*p(1)*( H - p(1) )/**sin( M_PI*FunctionTime::get_time()/8. )*//(H*H); return return_value; } - template inline Tensor<1,dim> Velocity::gradient( const Point &p, const unsigned int ) const{ - Tensor<1,dim> return_value; - switch( MultiComponentFunction::comp ){ - case 0: - return_value[0] = 0.; - return_value[1] = -std::cos( Function::get_time() ); - break; - case 1: - return_value[0] = std::cos( Function::get_time() ); - return_value[1] = 0.; - break; - default: - Assert( false, ExcNotImplemented() ); - }; - return return_value; - } - - template void Velocity::gradient_list( const std::vector > &points, - std::vector< Tensor<1,dim> > &gradients, - const unsigned int ) const{ - const unsigned int n_points = points.size(); - Assert( gradients.size() == n_points, ExcDimensionMismatch( gradients.size(), n_points ) ); - for( unsigned int i=0; i::gradient( points[i] ); - } - template class Pressure: public Function{ public: Pressure( const double initial_time = 0.0 ); virtual double value( const Point &p, const unsigned int component = 0 ) const; - virtual Tensor<1,dim> gradient( const Point &p, const unsigned int component = 0 ) const; virtual void value_list( const std::vector< Point > &points, std::vector &values, const unsigned int component = 0 ) const; - virtual void gradient_list( const std::vector< Point > &points, std::vector< Tensor<1,dim> > &gradients, - const unsigned int component = 0 ) const; }; template Pressure::Pressure( const double initial_time ): Function( 1, initial_time ){} template inline double Pressure::value( const Point &p, const unsigned int ) const{ - return std::sin( p(0) )*std::sin( p(1) )*std::sin( Function::get_time() ); - } - - template inline Tensor<1,dim> Pressure::gradient( const Point &p, const unsigned int ) const{ - return Point( std::cos( p(0) )*std::sin( p(1) )*std::sin( Function::get_time() ), - std::sin( p(0) )*std::cos( p(1) )*std::sin( Function::get_time() ) ); + return 0.; } template void Pressure::value_list( const std::vector > &points, @@ -364,16 +269,7 @@ namespace EquationData{ values[i] = Pressure::value( points[i] ); } - template inline void Pressure::gradient_list( const std::vector > &points, - std::vector< Tensor<1,dim> > &gradients, - const unsigned int ) const{ - const unsigned int n_points = points.size(); - Assert( gradients.size() == n_points, ExcDimensionMismatch( gradients.size(), n_points ) ); - for (unsigned int i=0; i::gradient( points[i] ); - } - - template class Force: public MultiComponentFunction{ +/* template class Force: public MultiComponentFunction{ public: Force( const double initial_time =0.0 ); virtual double value( const Point &p, const unsigned int component = 0 ) const; @@ -393,27 +289,9 @@ namespace EquationData{ } template inline double Force::value( const Point &p, const unsigned int ) const{ - double t = Function::get_time(), - cosx = std::cos( p(0) ), - sinx = std::sin( p(0) ), - cosy = std::cos( p(1) ), - siny = std::sin( p(1) ), - cost = std::cos(t), - sint = std::sin(t), - return_value = 0.; - switch( MultiComponentFunction::comp ){ - case 0: - return_value = p(1)*sint - p(0)*cost*cost + cosx*siny*sint; - break; - case 1: - return_value = -p(0)*sint - p(1)*cost*cost + sinx*cosy*sint ; - - break; - default: - Assert( false, ExcNotImplemented() ); - }; + double return_value = 0.; return return_value; - } + }*/ } @@ -427,11 +305,6 @@ template class Navier_Stokes_Projection{ Navier_Stokes_Projection( const RunTimeParameters::Data_Storage &data ); ~Navier_Stokes_Projection(); void run( const bool verbose = false, const unsigned int n_of_plots = 10 ); -/// - void Initialize(); - void set_dt( const double ddt ); - void Post_Process(); -/// protected: RunTimeParameters::Method_Formulation type; @@ -439,10 +312,11 @@ template class Navier_Stokes_Projection{ double dt; double t_0, T, Re; - EquationData::Force rhs; +// EquationData::Force rhs; EquationData::Velocity vel_exact; std::map boundary_values; + std::vector boundary_indicators; Triangulation triangulation; DoFHandler dof_handler_velocity, dof_handler_pressure; @@ -462,12 +336,11 @@ template class Navier_Stokes_Projection{ SparseILU prec_velocity[dim]; SparseDirectUMFPACK prec_mass, prec_pressure; - ConvergenceTable convergence_table; - DeclException2( ExcInvalidTimeStep, double, double, <<" The time step "< class Navier_Stokes_Projection{ private: unsigned int vel_max_its, vel_Krylov_size, vel_off_diagonals, vel_update_prec; double vel_eps, vel_diag_strength; - unsigned int proj_max_its, proj_off_diagonals; - double proj_eps, proj_diag_strength; - - unsigned int pres_max_its, pres_off_diagonals; - double pres_eps, pres_diag_strength; inline void init_velocity_matrices(); inline void init_pressure_matrices(); inline void init_gradient_operator(); typedef std_cxx1x::tuple< typename DoFHandler::active_cell_iterator, - typename DoFHandler::active_cell_iterator - > IteratorTuple; + typename DoFHandler::active_cell_iterator + > IteratorTuple; typedef parallel::internal::SynchronousIterators SIterators; struct InitGradPerTaskData{ unsigned int d, vel_dpc, pres_dpc; @@ -552,31 +420,23 @@ template Navier_Stokes_Projection::~Navier_Stokes_Projection(){ dof_handler_pressure.clear(); } -template void Navier_Stokes_Projection::set_dt( const double ddt ){ - AssertThrow( not ( ( ddt <= 0. ) or ( ddt > .5*T ) ), ExcInvalidTimeStep( ddt, .5*T ) ); - dt = ddt; -} // @sect4{ Navier_Stokes_Projection::Navier_Stokes_Projection } // In the constructor, we just read all the data from the Data_Storage -// object that is passed as an argument, verify that the read data is reasonable +// object that is passed as an argument, verify that the data we read is reasonable // and, finally, create the triangulation and load the initial data. template Navier_Stokes_Projection::Navier_Stokes_Projection( const RunTimeParameters::Data_Storage &data ): - type( data.form ), deg( data.pressure_degree ), dt( data.initial_dt ), t_0( data.initial_time ), - T( data.final_time ), Re( data.Reynolds ), rhs( data.initial_time ), + type( data.form ), deg( data.pressure_degree ), dt( data.dt ), t_0( data.initial_time ), + T( data.final_time ), Re( data.Reynolds ), /*rhs( data.initial_time ),*/ vel_exact( data.initial_time ), dof_handler_velocity( triangulation ), dof_handler_pressure( triangulation ), fe_velocity( deg+1 ), fe_pressure( deg ), quadrature_pressure( deg+1 ), quadrature_velocity( deg+2 ), vel_max_its( data.vel_max_iterations ), vel_Krylov_size( data.vel_Krylov_size ), vel_off_diagonals( data.vel_off_diagonals ), vel_update_prec( data.vel_update_prec ), vel_eps( data.vel_eps ), - vel_diag_strength( data.vel_diag_strength), - proj_max_its( data.proj_max_iterations ), proj_off_diagonals( data.proj_off_diagonals ), - proj_eps( data.proj_eps ), proj_diag_strength( data.proj_diag_strength ), - pres_max_its( data.pres_max_iterations), pres_off_diagonals( data.pres_off_diagonals ), - pres_eps( data.pres_eps ), pres_diag_strength( data.pres_diag_strength ) + vel_diag_strength( data.vel_diag_strength) { if(deg < 1) std::cout<<" WARNING: The chosen pair of finite element spaces is not stable."< Navier_Stokes_Projection::Navier_Stokes_Projection( // degrees of freedom and renumbers them, and initializes the matrices and vectors // that we will use. template void Navier_Stokes_Projection::Create_Triangulation( const unsigned int n_of_refines ){ - GridGenerator::hyper_ball( triangulation ); - static const HyperBallBoundary boundary; - triangulation.set_boundary( 0, boundary ); + GridIn grid_in; + grid_in.attach_triangulation( triangulation ); + std::ifstream file( "nsbench2.inp" ); + Assert( file, ExcFileNotOpen( "nsbench2.inp" ) ); + grid_in.read_ucd(file); + file.close(); + + std::cout<<" Number of refines = "< void Navier_Stokes_Projection::Create_Triangulation( cons init_pressure_matrices(); init_gradient_operator(); - pres_n.reinit( dof_handler_pressure.n_dofs() ); pres_n_minus_1.reinit( dof_handler_pressure.n_dofs() ); phi_n.reinit( dof_handler_pressure.n_dofs() ); @@ -705,7 +571,7 @@ template void Navier_Stokes_Projection::init_pressure_matrices(){ // For the gradient operator, we start by initializing the sparsity pattern and compressing it. // It is important to notice here that the gradient operator acts from the pressure space // into the velocity space, so we have to deal with two different finite element spaces. To keep -// the loops synchronized, we use the typedef's that we have defined before,namely +// the loops synchronized, we use the typedef's that we have defined before, namely // PairedIterators and SIterators. template void Navier_Stokes_Projection::init_gradient_operator(){ spar_pattern_pres_vel.reinit( dof_handler_velocity.n_dofs(), dof_handler_pressure.n_dofs(), @@ -773,10 +639,17 @@ template void Navier_Stokes_Projection::copy_gradient_local_to_glo // and so it is by default set to false template void Navier_Stokes_Projection::run( const bool verbose, const unsigned int n_of_plots ){ unsigned int n_steps = ( T - t_0 )/dt; - rhs.set_time( 2.*dt ); +// rhs.set_time( 2.*dt ); vel_exact.set_time( 2.*dt ); + plot_solution(1); + for( unsigned int n = 2; n<=n_steps; ++n ){ + if( n%n_of_plots == 0 ){ + if( verbose ) + std::cout<<" Plotting Solution"< void Navier_Stokes_Projection::interpolate_velocity(){ @@ -821,8 +690,8 @@ template void Navier_Stokes_Projection::diffusion_step( const bool assemble_advection_term(); for( unsigned int d=0; d void Navier_Stokes_Projection::diffusion_step( const bool vel_it_matrix[d].add( 1., vel_Advection ); vel_exact.set_component(d); - VectorTools::interpolate_boundary_values( dof_handler_velocity, 0, vel_exact, boundary_values ); + std::vector::const_iterator boundaries = boundary_indicators.begin(), + b_end = boundary_indicators.end(); + for( ; boundaries not_eq b_end; ++boundaries ){ + switch( *boundaries ){ + case 1: + VectorTools::interpolate_boundary_values( dof_handler_velocity, *boundaries, + ZeroFunction(), boundary_values ); + break; + case 2: + VectorTools::interpolate_boundary_values( dof_handler_velocity, *boundaries, + vel_exact, boundary_values ); + break; + /// + case 3: + if( d not_eq 0 ) + VectorTools::interpolate_boundary_values( dof_handler_velocity, *boundaries, + ZeroFunction(), boundary_values ); + break; + case 4: + VectorTools::interpolate_boundary_values( dof_handler_velocity, *boundaries, + ZeroFunction(), boundary_values ); + break; + default: + Assert( false, ExcNotImplemented() ); + } + } MatrixTools::apply_boundary_values( boundary_values, vel_it_matrix[d], u_n[d], force[d] ); } + Threads::TaskGroup tasks; for(unsigned int d=0; d::AdditionalData( vel_diag_strength, vel_off_diagonals ) ); - tasks += Threads::new_task( &Navier_Stokes_Projection::diffusion_component_solve, *this, d ); + tasks += Threads::new_task( &Navier_Stokes_Projection::diffusion_component_solve, *this, d ); } tasks.join_all(); } -template void Navier_Stokes_Projection::diffusion_component_solve( const unsigned int d){ +template void Navier_Stokes_Projection::diffusion_component_solve( const unsigned int d ){ SolverControl solver_control( vel_max_its, vel_eps*force[d].l2_norm() ); - SolverGMRES<> gmres( solver_control, SolverGMRES<>::AdditionalData() ); + SolverGMRES<> gmres( solver_control, SolverGMRES<>::AdditionalData( vel_Krylov_size ) ); gmres.solve( vel_it_matrix[d], u_n[d], force[d], prec_velocity[d] ); } @@ -906,7 +801,6 @@ template void Navier_Stokes_Projection::copy_advection_local_to_gl - // @sect4{Navier_Stokes_Projection::projection_step} // This implements the projection step. template void Navier_Stokes_Projection::projection_step( const bool reinit_prec ){ @@ -957,9 +851,8 @@ template void Navier_Stokes_Projection::update_pressure( const boo // @sect4{ Navier_Stokes_Projection::plot_solution } -// At this stage, we only output the vorticity of the flow. This only works in 2d and -// WILL be changed. -/// +// This method plots the current solution. It is an adaptation of +// step-31 **** and so I will not elaborate on it. template void Navier_Stokes_Projection::plot_solution( const unsigned int step ){ const FESystem joint_fe( fe_velocity, dim, fe_pressure, 1 ); DoFHandler joint_dof_handler( triangulation ); @@ -982,15 +875,12 @@ template void Navier_Stokes_Projection::plot_solution( const unsig for( unsigned int i=0; i void Navier_Stokes_Projection::plot_solution( const unsig Assert( false, ExcInternalError() ); } } - std::vector joint_solution_names( dim, "v" ); joint_solution_names.push_back( "p" ); - DataOut data_out; data_out.attach_dof_handler (joint_dof_handler); - std::vector< DataComponentInterpretation::DataComponentInterpretation > component_interpretation( dim+1, DataComponentInterpretation::component_is_part_of_vector ); component_interpretation[dim] = DataComponentInterpretation::component_is_scalar; - data_out.add_data_vector( joint_solution, joint_solution_names, DataOut::type_dof_data, component_interpretation ); - data_out.build_patches( deg + 1 ); - std::ostringstream filename; filename<<"solution-"<Navier_Stokes_Projection::Post_Process} -// Having reached the final time T, we want to measure the error that we have made. -// This method is responsible for that. Saves the results in a ConvergenceTable -// object which later we can print or compute things with it.
-// The way we compute the errors is very similar to previous tutorials. However, we need the -// pressure to have mean value zero, so we compute its mean value and subtract it from the computed -// pressure. -template void Navier_Stokes_Projection::Post_Process(){ - double tmp, vel_err_L2=0., vel_err_H1=0., pres_err_L2; - - Vector differences( triangulation.n_active_cells() ); - - vel_exact.set_time(T); - for( unsigned int d=0; d pres_exact(T); - differences = 0.; - VectorTools::integrate_difference( dof_handler_pressure, pres_n, pres_exact, - differences, quadrature_pressure, VectorTools::L2_norm ); - pres_err_L2 = differences.l2_norm(); - - convergence_table.add_value( "dt" , dt ); - convergence_table.add_value( "u_L2" , vel_err_L2 ); - convergence_table.add_value( "u_H1" , vel_err_H1 ); - convergence_table.add_value( "pres_L2", pres_err_L2 ); - - convergence_table.set_precision( "dt" , 5 ); - convergence_table.set_precision( "u_L2" , 5 ); - convergence_table.set_precision( "u_H1" , 5 ); - convergence_table.set_precision( "pres_L2", 5 ); - - convergence_table.set_scientific( "u_L2" , true ); - convergence_table.set_scientific( "u_H1" , true ); - convergence_table.set_scientific( "pres_L2", true ); - - convergence_table.write_text(std::cout); -} - - // @sect3{ The main function } // The main function looks very much like in all the other tutorial programs. @@ -1089,14 +915,7 @@ int main(){ data.read_data( "parameter-file.prm" ); deallog.depth_console( data.verbose?2:0 ); Navier_Stokes_Projection<2> test( data ); - for( double dt = data.initial_dt; dt >= data.final_dt; dt /= data.dt_decrement ){ - std::cout<<" dt = "<