+++ /dev/null
-/*
- Header file for the classes that define the exact solution
- and driving force.
-
- These classes are much like any other Function class
- so we are not going to comment on them
-
- by Abner Salgado.
-*/
-#ifndef _EQ_DATA_H_
-#define _EQ_DATA_H_
-
-
-/*
- We want to be able to select which component of vector
- valued functions we are going to work on
-*/
-#include "../include/MultiComponentFunction.h"
-
-
-// this is dealii
-#include <base/point.h>
-
-
-// We need to define sines and cosines
-#include <cmath>
-
-
-
-// This is ugly, we have to remove it
-const double PI = std::acos( -1. );
-
-
-
-// The Velocity function
-template<int dim> class Velocity: public Multi_Component_Function<dim> {
- public:
- Velocity(const double initial_time =0.0);
- virtual double value(const Point<dim> &p, const unsigned int component = 0) const;
- virtual Tensor<1,dim> gradient(const Point<dim> &p, const unsigned int component=0) const;
- virtual void value_list( const std::vector< Point<dim> > &points, std::vector<double> &values,
- const unsigned int component = 0 ) const;
- virtual void gradient_list( const std::vector< Point<dim> > &points, std::vector< Tensor<1,dim> > &gradients,
- const unsigned int component = 0 ) const;
-};
-
-
-
-// The Pressure function
-template<int dim> class Pressure: public Function<dim>{
- public:
- Pressure(const double initial_time = 0.0);
- virtual double value(const Point<dim> &p, const unsigned int component = 0) const;
- virtual Tensor<1,dim> gradient(const Point<dim> &p, const unsigned int component=0) const;
- virtual void value_list( const std::vector< Point<dim> > &points, std::vector<double> &values,
- const unsigned int component = 0 ) const;
- virtual void gradient_list( const std::vector< Point<dim> > &points, std::vector< Tensor<1,dim> > &gradients,
- const unsigned int component = 0 ) const;
-};
-
-
-
-// The Force function
-template<int dim> class Force: public Multi_Component_Function<dim>{
- public:
- Force( const double initial_time =0.0 );
- virtual double value( const Point<dim> &p, const unsigned int component = 0 ) const;
- virtual void value_list( const std::vector< Point<dim> > &points, std::vector<double> &values, const unsigned int component = 0 ) const;
-};
-
-#endif
+++ /dev/null
-// This class is supposed to aid in the reading of parameter files
-#ifndef _FILE_READER_H_
-#define _FILE_READER_H_
-
-
-
-#include <base/parameter_handler.h>
-
-
-
-#include <fstream>
-#include <iostream>
-
-
-
-using namespace dealii;
-
-enum Method_Formulation{
- METHOD_STANDARD,
- METHOD_ROTATIONAL
-};
-
-class Data_Storage{
- public:
- Data_Storage();
- ~Data_Storage();
- void read_data( char *filename );
- void print_usage();
- void print_status() const;
- // The data itself
- //// The type of method we want to use
- Method_Formulation form;
- //// physical data
- double initial_time,
- final_time,
- Reynolds;
- //// Time stepping data
- double initial_dt,
- final_dt,
- dt_decrement;
- //// Space discretization data
- unsigned int n_of_global_refines,
- pressure_degree;
- //// Data to solve the velocity
- unsigned int vel_max_iterations,
- vel_Krylov_size,
- vel_off_diagonals,
- vel_update_prec;
- double vel_eps,
- vel_diag_strength;
- //// Data to solve the projection
- unsigned int proj_max_iterations,
- proj_off_diagonals;
- double proj_eps,
- proj_diag_strength;
- //// Data to do the pressure update step
- unsigned int pres_max_iterations,
- pres_off_diagonals;
- double pres_eps,
- pres_diag_strength;
- //// Verbosity
- bool verbose;
- //// Frequency of the outputted data
- unsigned int output;
-
- protected:
- ParameterHandler prm;
-};
-
-#endif
+++ /dev/null
-/*
- Header file for Multi Component Functions.
- The whole purpose of this class to exist is so we do not have to write the function
- void set_component( const unsigned int d )
- twice
-
- by Abner Salgado.
-*/
-#ifndef _MULTI_COMPONENT_FUNCTION_H_
-#define _MULTI_COMPONENT_FUNCTION_H_
-
-
-/*
- This is basically a workaround for vector valued functions
- when we just want one of its components
- so they are derived from the Function class
-*/
-#include <base/function.h>
-
-
-using namespace dealii;
-
-
-/*
- The only funcitonality this class has is that it provides a common wrapper
- for vector valued functions for the case when you only want to ask them for
- one of their components
-*/
-template<int dim> class Multi_Component_Function: public Function<dim> {
- public:
- /*
- Constructor. It does not do anything interesting but set the initial
- time for the function to evaluate its values
- */
- Multi_Component_Function( const double initial_time = 0. );
- /*
- This is the whole reason this class exists.
- It wouldn't seem logical to have this written twice. One for the
- Velocity function and one for the Force function.
- So we are instead going to derive these from this class
- and so adding this functionality of selecting which component you
- want to give the value of
- */
- void set_component(const unsigned int d );
- protected:
- // The current component we are working on
- unsigned int component;
-};
-
-
-#endif
+++ /dev/null
-/*
- Definition of the Projection Method for constant density
- Navier Stokes
- by Abner Salgado.
-*/
-#ifndef _NAVIER_STOKES_H_
-#define _NAVIER_STOKES_H_
-
-
-
-
-/*
- First we include the local files, where we have the definitions of
- 1.- The exact solution
- 2.- The class that reads runtime parameters from a file
-*/
-#include "EqData.h"
-#include "FileReader.h"
-
-
-
-/*
- These includes are to get the quadrature
- and to handle the convergence table
-*/
-#include <base/quadrature_lib.h>
-#include <base/convergence_table.h>
-
-
-
-/*
- The includes needed to manage multithreading
-*/
-#include <base/multithread_info.h>
-#include <base/thread_management.h>
-
-
-
-/*
- Grid management includes.
- Triangulation handler
- Grid generator & refinement, accessor and iterators
-*/
-#include <grid/tria.h>
-#include <grid/grid_generator.h>
-#include <grid/grid_refinement.h>
-#include <grid/tria_accessor.h>
-#include <grid/tria_iterator.h>
-#include <grid/tria_boundary_lib.h>
-#include <grid/grid_in.h>
-
-
-
-/*
- Degrees of freedom management includes
- Handling, accessing, computing useful stuff from and for them and renumbering.
-*/
-#include <dofs/dof_handler.h>
-#include <dofs/dof_accessor.h>
-#include <dofs/dof_tools.h>
-#include <dofs/dof_renumbering.h>
-#include <dofs/dof_constraints.h>
-
-
-
-/*
- Finite element management includes
- Definition of finite elements and extraction of info from them
- and from finite element functions.
-*/
-#include <fe/fe_q.h>
-#include <fe/fe_values.h>
-#include <fe/fe_tools.h>
-
-
-
-/*
- Linear algebra includes
- Sparse matrices, vectors, preconditioners and solvers.
-*/
-#include <lac/vector.h>
-#include <lac/sparse_matrix.h>
-#include <lac/solver_cg.h>
-#include <lac/precondition.h>
-#include <lac/solver_gmres.h>
-#include <lac/sparse_ilu.h>
-#include <lac/sparse_direct.h>
-
-
-
-/*
- Numerical algorithms includes
- Assembly of right hand sides, computation of errors, output of the solution, etc
-*/
-#include <numerics/matrices.h>
-#include <numerics/vectors.h>
-#include <numerics/data_out.h>
-
-
-/*
- This class is the core of all the program
- It is the implementation of the projection method for the Navier-Stokes equations.
-*/
-template<int dim> class Navier_Stokes_Projection{
- public:
- /*
- Constructor. It takes as a parameter a reference to a Data_Storage class, which is
- basically a container of all the data this class needs to be defined (and some other stuff)
- Here we basically just copy the needed data and/or assign the values that are needed.
- */
- Navier_Stokes_Projection( const Data_Storage &data );
- /*
- Destructor. Completely trivial.
- */
- ~Navier_Stokes_Projection();
- /*
- This function creates the triangulation and refines it the required number of times
- It also initializes all the quantities that are mesh dependent but not time stepping dependent
- i.e. the size of matrices and vectors.
- */
- void Create_Triangulation( const unsigned int n_of_refines );
- /*
- Having created a mesh and once the time step is set, the next step is to initialize the method
- i.e. compute the matrices that are never going to change, load the initial data, etc.
- This is what this function does.
- */
- void Initialize();
- /*
- This is the time marching function, which starting at t_0 advances in time using the projection method
- with time step dt until T.
- The boolean parameter that it takes is to enable information about what the function is doing at the present
- moment, i.e. diffusion, projection substep; updating preconditioners etc.
- This is useful mostly for debuggin purposes and so it is by default set to false
- */
- void run( const bool verbose = false, const unsigned int n_of_plots = 10 );
- /*
- 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
- */
- void Post_Process();
- /*
- The whole reason for this class, and all this code after all, is to make convergence tests
- with respect to the time discretization. For this reason we need to be able to vary the time step
- we are going to work with. This method sets up the time step that is to be used in the given run
- */
- void set_dt( const double ddt );
-
- protected:
- // The type of method
- Method_Formulation type;
-
- // Discretization data
- //// The polynomial degree
- unsigned int deg;
- //// The time step
- double dt;
-
- // Physical data:
- //// initial time, final time, Reynolds number
- double t_0, T, Re;
- //// The external driving force
- Force<dim> rhs;
- //// The exact velocity (to apply boundary values)
- Velocity<dim> vel_exact;
- //// The boundary conditions
- std::map<unsigned int, double> boundary_values;
-
- // Finite Element Spaces data
- //// The mesh
- Triangulation<dim> triangulation;
- //// The DoF handlers
- DoFHandler<dim> dof_handler_velocity, dof_handler_pressure;
- //// The polynomial spaces
- FE_Q<dim> fe_velocity, fe_pressure;
- //// The quadrature formulae
- QGauss<dim> quadrature_pressure, quadrature_velocity;
-
- /*
- Linear Algebra Data
- The sparsity patterns where the matrices will live
- */
- SparsityPattern spar_pattern_velocity, spar_pattern_pressure, spar_pattern_pres_vel;
- /*
- The actual matrices. The projection matrix never changes.
- For the velocity a part of the matrix never changes (neither in time nor for each component),
- namely the part related with the time derivative and the diffusion
- but, the advection part changes and so we need dim+1 matrices. One to store
- the constant part, and dim for each matrix at each iteration
- */
- SparseMatrix<double> vel_Laplace_plus_Mass, vel_it_matrix[dim], vel_Mass, vel_Laplace,
- pres_Laplace, pres_Mass, pres_Diff[dim];
- /*
- We need to regularize the Laplace operator for the pressure
- for that we use a constraint
- */
- ConstraintMatrix pres_regularization;
- //// The solutions at times n and (n-1)
- Vector<double> pres_n, pres_n_minus_1, phi_n, phi_n_minus_1, u_n[dim], u_n_minus_1[dim],
- //// The time extrapolation of the velocity, used to make the semi-implicit time discretization
- //// of the advection term
- u_star[dim],
- //// Right hand side for the component of the momentum equation
- force[dim],
- //// Temporary arrays
- v_tmp, pres_tmp;
- //// The preconditioners
-// SparseILU<double> prec_velocity[dim];
- SparseDirectUMFPACK prec_velocity[dim], prec_mass, prec_pressure;
-
- // The convergence table
- ConvergenceTable convergence_table;
-
- // Exception we will throw if the time step is invalid
- DeclException2( ExcInvalidTimeStep, double, double, <<" The time step "<<arg1<<" is out of range."<<std::endl
- <<" The permitted range is (0,"<<arg2<<"]");
-
- /*
- A diffusion step. The boolean parameter that it takes indicates if we want the preconditioners
- for each component of the velocity to be updated. Recall that the operator that these matrices
- represent is of the form
- 1.5/dt Mass + 1/Re Laplace_h + ( u^*.D )
- so the dominating part is
- 1.5/dt Mass + 1/Re Laplace_h
- and we can recycle the preconditioner for several steps to save computational time
- */
- inline void diffusion_step( const bool reinit_prec );
- /*
- A projection step. The parameter that it takes indicates if we want the preconditioner to be
- updated. Since this matrix never changes, ideally this should be done just once.
- */
- inline void projection_step( const bool reinit_prec );
- /*
- A pressure update step.
- Capable of recognizing which formulation of the method we are using.
- If it is the standard formulation, then the pressure update is done by
- p^{n+1} = p^{n} + phi^{n+1}
- which involves just adding numerical values.
- On the other hand, for the rotational formulation of the method the pressure update
- is done by the formula
- p^{n+1} = p^{n} + phi^{n+1} +1/Re div u^{n+1}
- which in any case involves the solution of a mass matrix problem.
- NOTE: Here we assume that the mesh is uniformly refined, and thus we do not use
- a preconditioner for this mass matrix problem
- */
- inline void update_pressure( const bool reinit_prec );
- private:
-
- // A bunch of data for the solvers and preconditioners
- // data solve velocity
- unsigned int vel_max_its, vel_Krylov_size, vel_off_diagonals, vel_update_prec;
- double vel_eps, vel_diag_strength;
- // data solve projection
- unsigned int proj_max_its, proj_off_diagonals;
- double proj_eps, proj_diag_strength;
- // data solve pressure update
- unsigned int pres_max_its, pres_off_diagonals;
- double pres_eps, pres_diag_strength;
-
- // Extremely simple function that computes the velocity extrapolant. Added for readability
- inline void interpolate_velocity();
-
- //////// initialization of constant wrt dt matrices
- inline void init_velocity_matrices();
- inline void init_pressure_matrices();
- inline void init_gradient_operator();
-
- // Assembly of the advection term
- inline void assemble_advection_term( const unsigned int d );
-
- // Solution of each component of the velocity on the diffusion step
- inline void threaded_solve( const unsigned int d );
- inline void plot_solution( const unsigned int step );
-};
-
-
-
-#endif
+++ /dev/null
-/*
- Implementation of the classes that represent the exact solution
- and external force FOR THE PROBLEM ON THE DISK
-
-
- There's not much to say here. All these classes just represent
- mathematical formulae
-
- by Abner Salgado.
-*/
-
-
-
-#include "../include/EqData.h"
-
-
-
-// Force function methods
-template<int dim> Force<dim>::Force( const double initial_time ): Multi_Component_Function<dim>( initial_time ){
-}
-
-
-
-template<int dim> void Force<dim>::value_list( const std::vector<Point<dim> > &points, std::vector<double> &values,
- const unsigned int ) const{
- const unsigned int n_points = points.size();
- Assert( values.size() == n_points, ExcDimensionMismatch( values.size(), n_points ) );
- for (unsigned int i=0; i<n_points; ++i)
- values[i] = Force<dim>::value( points[i] );
-}
-
-
-
-template<int dim> inline double Force<dim>::value(const Point<dim> &p, const unsigned int) const{
- double t = FunctionTime::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( Multi_Component_Function<dim>::component ){
- case 0:
- // y*sin(t)-x*cos(t)^2+cos(x)*sin(y)*sin(t)
- return_value = p(1)*sint - p(0)*cost*cost + cosx*siny*sint;
- break;
- case 1:
- // -x*sin(t)-y*cos(t)^2+sin(x)*cos(y)*sin(t)
- return_value = -p(0)*sint - p(1)*cost*cost + sinx*cosy*sint ;
-
- break;
- default:
- Assert( false, ExcNotImplemented() );
- };
- return return_value;
-}
-
-
-
-// Velocity function methods
-template<int dim> Velocity<dim>::Velocity(const double initial_time): Multi_Component_Function<dim>( initial_time ){
-}
-
-
-
-template<int dim> void Velocity<dim>::value_list( const std::vector<Point<dim> > &points, std::vector<double> &values,
- const unsigned int ) const{
- const unsigned int n_points = points.size();
- Assert( values.size() == n_points, ExcDimensionMismatch( values.size(), n_points ) );
- for (unsigned int i=0; i<n_points; ++i)
- values[i] = Velocity<dim>::value( points[i] );
-}
-
-
-
-template<int dim> inline double Velocity<dim>::value(const Point<dim> &p, const unsigned int) const{
- double return_value = std::cos( Function<dim>::get_time() );
- switch( Multi_Component_Function<dim>::component ){
- case 0:
- // -y*cos(t)
- return_value *= -p(1);
- break;
- case 1:
- // x*cos(t)
- return_value *= p(0);
- break;
- default:
- Assert( false, ExcNotImplemented() );
- };
- return return_value;
-}
-
-
-
-template<int dim> inline Tensor<1,dim> Velocity<dim>::gradient(const Point<dim> &p, const unsigned int) const{
- Tensor<1,dim> return_value;
- switch( Multi_Component_Function<dim>::component ){
- // [0, -cos(t)]
- case 0:
- return_value[0] = 0.;
- return_value[1] = -std::cos( Function<dim>::get_time() );
- break;
- case 1:
- // [cos(t), 0]
- return_value[0] = std::cos( Function<dim>::get_time() );
- return_value[1] = 0.;
- break;
- default:
- Assert( false, ExcNotImplemented() );
- };
- return return_value;
-}
-
-
-
-template<int dim> void Velocity<dim>::gradient_list( const std::vector<Point<dim> > &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<n_points; ++i)
- gradients[i] = Velocity<dim>::gradient( points[i] );
-}
-
-
-
-// Pressure function methods
-template<int dim> Pressure<dim>::Pressure(const double initial_time): Function<dim>(1,initial_time){}
-
-
-
-template<int dim> inline double Pressure<dim>::value(const Point<dim> &p, const unsigned int) const{
- // sin(x-y+t)
- return std::sin( p(0) )*std::sin( p(1) )*std::sin( Function<dim>::get_time() );
-}
-
-
-
-template<int dim> inline Tensor<1,dim> Pressure<dim>::gradient(const Point<dim> &p, const unsigned int) const{
- // [cos(x)*sin(y)*sin(t), sin(x)*cos(y)*sin(t)]
- return Point<dim>( std::cos( p(0) )*std::sin( p(1) )*std::sin( Function<dim>::get_time() ),
- std::sin( p(0) )*std::cos( p(1) )*std::sin( Function<dim>::get_time() ) );
-}
-
-
-
-template<int dim> void Pressure<dim>::value_list( const std::vector<Point<dim> > &points, std::vector<double> &values,
- const unsigned int ) const{
- const unsigned int n_points = points.size();
- Assert( values.size() == n_points, ExcDimensionMismatch( values.size(), n_points ) );
- for (unsigned int i=0; i<n_points; ++i)
- values[i] = Pressure<dim>::value( points[i] );
-}
-
-
-
-template<int dim> inline void Pressure<dim>::gradient_list( const std::vector<Point<dim> > &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<n_points; ++i)
- gradients[i] = Pressure<dim>::gradient( points[i] );
-}
-
-
-
-// explicit template instantiation
-template class Force<deal_II_dimension>;
-template class Velocity<deal_II_dimension>;
-template class Pressure<deal_II_dimension>;
+++ /dev/null
-/*
- Implementation of the classes that represent the exact solution
- and external force.
-
- There's not much to say here. All these classes just represent
- mathematical formulae
-
- by Abner Salgado.
-*/
-
-
-
-#include "../include/EqData.h"
-
-
-
-// Force function methods
-template<int dim> Force<dim>::Force( const double initial_time ): Multi_Component_Function<dim>( initial_time ){
-}
-
-
-
-template<int dim> void Force<dim>::value_list( const std::vector<Point<dim> > &points, std::vector<double> &values,
- const unsigned int ) const{
- const unsigned int n_points = points.size();
- Assert( values.size() == n_points, ExcDimensionMismatch( values.size(), n_points ) );
- for (unsigned int i=0; i<n_points; ++i)
- values[i] = Force<dim>::value( points[i] );
-}
-
-
-
-template<int dim> inline double Force<dim>::value(const Point<dim> &p, const unsigned int) const{
- double t = FunctionTime::get_time(),
- sin_2pi_x = std::sin( 2.*PI*p(0) ),
- sin_2pi_y = std::sin( 2.*PI*p(1) ),
- cos_2pi_x = std::cos( 2.*PI*p(0) ),
- cos_2pi_y = std::cos( 2.*PI*p(1) ),
- sin_pi_x = std::sin( PI*p(0) ),
- sin_pi_y = std::sin( PI*p(1) ),
- cos_pi_x = std::cos( PI*p(0) ),
- cos_pi_y = std::cos( PI*p(1) ),
- cos_t = std::cos( t ),
- sin_t = std::sin( t ),
- return_value = 0.;
- switch( Multi_Component_Function<dim>::component ){
- case 0:
- // Pi*sin(2*Pi*y)*sin(Pi*x)^2*cos(t)
- // -2*Pi^3*sin(2*Pi*y)*cos(Pi*x)^2*sin(t)
- // +6*Pi^3*sin(2*Pi*y)*sin(Pi*x)^2*sin(t)
- // +2*Pi^3*sin(2*Pi*y)^2*sin(Pi*x)^3*sin(t)^2*cos(Pi*x)
- // -2*Pi^3*sin(2*Pi*x)*sin(Pi*y)^2*sin(t)^2*cos(2*Pi*y)*sin(Pi*x)^2
- // -Pi*sin(Pi*x)*sin(Pi*y)*sin(t)
- return_value = PI*sin_2pi_y*sin_pi_x*sin_pi_x*cos_t
- - 2.*PI*PI*PI*sin_2pi_y*cos_pi_x*cos_pi_x*sin_t
- + 6.*PI*PI*PI*sin_2pi_y*sin_pi_x*sin_pi_x*sin_t
- + 2.*PI*PI*PI*sin_2pi_y*sin_2pi_y*sin_pi_x*sin_pi_x*sin_pi_x*sin_t*sin_t*cos_pi_x
- - 2.*PI*PI*PI*sin_2pi_x*sin_pi_y*sin_pi_y*sin_t*sin_t*cos_2pi_y*sin_pi_x*sin_pi_x
- - PI*sin_pi_x*sin_pi_y*sin_t;
-
- break;
- case 1:
- // -Pi*sin(2*Pi*x)*sin(Pi*y)^2*cos(t)
- // -6*Pi^3*sin(2*Pi*x)*sin(Pi*y)^2*sin(t)
- // +2*Pi^3*sin(2*Pi*x)*cos(Pi*y)^2*sin(t)
- // +Pi*cos(Pi*x)*cos(Pi*y)*sin(t)
- // -2*Pi^3*sin(2*Pi*y)*sin(Pi*x)^2*sin(t)^2*cos(2*Pi*x)*sin(Pi*y)^2
- // +2*Pi^3*sin(2*Pi*x)^2*sin(Pi*y)^3*sin(t)^2*cos(Pi*y)
- return_value = -PI*sin_2pi_x*sin_pi_y*sin_pi_y*cos_t
- - 6.*PI*PI*PI*sin_2pi_x*sin_pi_y*sin_pi_y*sin_t
- + 2.*PI*PI*PI*sin_2pi_x*cos_pi_y*cos_pi_y*sin_t
- + PI*cos_pi_x*cos_pi_y*sin_t
- - 2.*PI*PI*PI*sin_2pi_y*sin_pi_x*sin_pi_x*sin_t*sin_t*cos_2pi_x*sin_pi_y*sin_pi_y
- + 2.*PI*PI*PI*sin_2pi_x*sin_2pi_x*sin_pi_y*sin_pi_y*sin_pi_y*sin_t*sin_t*cos_pi_y;
-
- break;
- default:
- Assert( false, ExcNotImplemented() );
- };
- return return_value;
-}
-
-
-
-// Velocity function methods
-template<int dim> Velocity<dim>::Velocity(const double initial_time): Multi_Component_Function<dim>( initial_time ){
-}
-
-
-
-template<int dim> void Velocity<dim>::value_list( const std::vector<Point<dim> > &points, std::vector<double> &values,
- const unsigned int ) const{
- const unsigned int n_points = points.size();
- Assert( values.size() == n_points, ExcDimensionMismatch( values.size(), n_points ) );
- for (unsigned int i=0; i<n_points; ++i)
- values[i] = Velocity<dim>::value( points[i] );
-}
-
-
-
-template<int dim> inline double Velocity<dim>::value(const Point<dim> &p, const unsigned int) const{
- double sin_2pi_x = std::sin( 2.*PI*p(0) ),
- sin_2pi_y = std::sin( 2.*PI*p(1) ),
- sin_pi_x = std::sin( PI*p(0) ),
- sin_pi_y = std::sin( PI*p(1) ),
- sin_t = std::sin( Function<dim>::get_time() ),
- return_value = 0.;
- switch( Multi_Component_Function<dim>::component ){
- case 0:
- // pi*sin(2*pi*y)*sin(pi*x)^2*sin(t)
- return_value = PI*sin_2pi_y*sin_pi_x*sin_pi_x*sin_t;
- break;
- case 1:
- // -pi*sin(2*pi*x)*sin(pi*y)^2*sin(t)
- return_value = -PI*sin_2pi_x*sin_pi_y*sin_pi_y*sin_t;
- break;
- default:
- Assert( false, ExcNotImplemented() );
- };
- return return_value;
-}
-
-
-
-template<int dim> inline Tensor<1,dim> Velocity<dim>::gradient(const Point<dim> &p, const unsigned int) const{
- Tensor<1,dim> return_value;
- double sin_2pi_x = std::sin( 2.*PI*p(0) ),
- sin_2pi_y = std::sin( 2.*PI*p(1) ),
- cos_2pi_x = std::cos( 2.*PI*p(0) ),
- cos_2pi_y = std::cos( 2.*PI*p(1) ),
- sin_pi_x = std::sin( PI*p(0) ),
- sin_pi_y = std::sin( PI*p(1) ),
- cos_pi_x = std::cos( PI*p(0) ),
- cos_pi_y = std::cos( PI*p(1) ),
- sin_t = std::sin( Function<dim>::get_time() );
- switch( Multi_Component_Function<dim>::component ){
- // [2*Pi^2*sin(2*Pi*y)*sin(Pi*x)*sin(t)*cos(Pi*x),
- // 2*Pi^2*cos(2*Pi*y)*sin(Pi*x)^2*sin(t)]
- case 0:
- return_value[0] = 2.*PI*PI*sin_2pi_y*sin_pi_x*sin_t*cos_pi_x;
- return_value[1] = 2.*PI*PI*cos_2pi_y*sin_pi_x*sin_pi_x*sin_t;
- break;
- case 1:
- // [-2*Pi^2*cos(2*Pi*x)*sin(Pi*y)^2*sin(t),
- // -2*Pi^2*sin(2*Pi*x)*sin(Pi*y)*sin(t)*cos(Pi*y)]
- return_value[0] = -2.*PI*PI*cos_2pi_x*sin_pi_y*sin_pi_y*sin_t;
- return_value[1] = -2.*PI*PI*sin_2pi_x*sin_pi_y*sin_t*cos_pi_y;
- break;
- default:
- Assert( false, ExcNotImplemented() );
- };
- return return_value;
-}
-
-
-
-template<int dim> void Velocity<dim>::gradient_list( const std::vector<Point<dim> > &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<n_points; ++i)
- gradients[i] = Velocity<dim>::gradient( points[i] );
-}
-
-
-
-// Pressure function methods
-template<int dim> Pressure<dim>::Pressure(const double initial_time): Function<dim>(1,initial_time){}
-
-
-
-template<int dim> inline double Pressure<dim>::value(const Point<dim> &p, const unsigned int) const{
- // cos(pi*x)*sin(pi*y)*sin(t)
- return std::cos( PI*p(0) )*std::sin( PI*p(1) )*std::sin( Function<dim>::get_time() );
-}
-
-
-
-template<int dim> inline Tensor<1,dim> Pressure<dim>::gradient(const Point<dim> &p, const unsigned int) const{
- // || -sin(Pi*x)*Pi*sin(Pi*y)*sin(t) ||
- // || cos(Pi*x)*cos(Pi*y)*Pi*sin(t) ||
- return ( PI*std::sin( Function<dim>::get_time() ) )* Point<dim>( -std::sin( PI*p(0) )*std::sin( PI*p(1) ),
- std::cos( PI*p(0) )*std::cos( PI*p(1) ) ) ;
-}
-
-
-
-template<int dim> void Pressure<dim>::value_list( const std::vector<Point<dim> > &points, std::vector<double> &values,
- const unsigned int ) const{
- const unsigned int n_points = points.size();
- Assert( values.size() == n_points, ExcDimensionMismatch( values.size(), n_points ) );
- for (unsigned int i=0; i<n_points; ++i)
- values[i] = Pressure<dim>::value( points[i] );
-}
-
-
-
-template<int dim> inline void Pressure<dim>::gradient_list( const std::vector<Point<dim> > &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<n_points; ++i)
- gradients[i] = Pressure<dim>::gradient( points[i] );
-}
-
-
-
-// explicit template instantiation
-template class Force<deal_II_dimension>;
-template class Velocity<deal_II_dimension>;
-template class Pressure<deal_II_dimension>;
+++ /dev/null
-#include "../include/FileReader.h"
-
-
-
-// Constructor. Here we set up the
-// Data Format we are going to use
-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. " );
-
- // Physical data of the problem
- prm.enter_subsection( "Physical data" );
- prm.declare_entry( "initial_time", "0.", Patterns::Double( 0. ), " The initial time of the simulation. " );
- prm.declare_entry( "final_time", "1.", Patterns::Double( 0. ), " The final time of the simulation. " );
- prm.declare_entry( "Reynolds", "1.", Patterns::Double( 0. ), " The Reynolds number. " );
- prm.leave_subsection();
-
- // Time stepping data of the problem
- 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.leave_subsection();
-
- // Space discretization data
- prm.enter_subsection( "Space discretization" );
- prm.declare_entry( "n_of_refines", "5", Patterns::Integer( 1, 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. " );
- prm.leave_subsection();
-
- // Velocity solution data
- prm.enter_subsection( "Data solve velocity" );
- prm.declare_entry( "max_iterations", "1000", Patterns::Integer( 1, 1000 ), " The maximal number of iterations GMRES must make. " );
- prm.declare_entry( "eps", "1e-12", Patterns::Double( 0. ), " The stopping criterion. " );
- prm.declare_entry( "Krylov_size", "30", Patterns::Integer(1), " The size of the Krylov subspace to be used. " );
- prm.declare_entry( "off_diagonals", "60", Patterns::Integer(0), " The number of off-diagonal elements ILU must compute. " );
- prm.declare_entry( "diag_strength", "0.01", Patterns::Double( 0. ), " Diagonal strengthening coefficient. " );
- prm.declare_entry( "update_prec", "15", Patterns::Integer(1), " This number indicates how often we need to update the preconditioner" );
- prm.leave_subsection();
-
- // Projection step data
- 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();
-
- // Pressure update data
- 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();
-
- // Verbosity of output
- prm.declare_entry( "verbose", "true", Patterns::Bool(), " This indicates whether the output of the solution process should be verbose. " );
-
- // How often we want the data to be outputted
- prm.declare_entry( "output", "10", Patterns::Integer(1), " This indicates between how many time steps we print the solution. " );
-}
-
-
-
-// Destructor. Does nothing
-Data_Storage::~Data_Storage(){
-}
-
-
-
-// Here is where all happens. We read all the data from the indicated file
-void Data_Storage::read_data( char *filename ){
- std::ifstream file;
- file.open( filename );
- if( not file )
- throw ExcFileNotOpen( filename );
- prm.read_input( file );
-
- std::string token = prm.get( "Method_Form" );
- if( token == std::string("rotational") )
- form = METHOD_ROTATIONAL;
- else
- form = METHOD_STANDARD;
-
- // Physical data of the problem
- prm.enter_subsection( "Physical data" );
- initial_time = prm.get_double( "initial_time" );
- final_time = prm.get_double( "final_time" );
- Reynolds = prm.get_double( "Reynolds" );
- prm.leave_subsection();
-
- // Time stepping data of the problem
- 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" );
- prm.leave_subsection();
-
- // Space discretization data
- prm.enter_subsection( "Space discretization" );
- n_of_global_refines = prm.get_integer( "n_of_refines" );
- pressure_degree = prm.get_integer( "pressure_fe_degree" );
- prm.leave_subsection();
-
- // Velocity solution data
- prm.enter_subsection( "Data solve velocity" );
- vel_max_iterations = prm.get_double( "max_iterations" );
- vel_eps = prm.get_double( "eps" );
- vel_Krylov_size = prm.get_integer( "Krylov_size" );
- vel_off_diagonals = prm.get_integer( "off_diagonals" );
- vel_diag_strength = prm.get_double( "diag_strength" );
- vel_update_prec = prm.get_integer( "update_prec" );
- prm.leave_subsection();
-
- // Projection step data
- 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();
-
- // Pressure update data
- 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();
-
- // Verbosity
- verbose = prm.get_bool( "verbose" );
-
- // Output frequency
- output = prm.get_integer( "output" );
-
- file.close();
-}
-
-
-
-// Prints the current values of the data
-// Mostly (if not only) used for debugging purposes
-void Data_Storage::print_status() const{
- std::cout<<"Method Form = "<<( (form == METHOD_ROTATIONAL)?"rotational":"standard" )<<std::endl
- <<"Physical Data:"<<std::endl
- <<" initial_time = "<<initial_time<<std::endl
- <<" final_time = "<<final_time<<std::endl
- <<" Re = "<<Reynolds<<std::endl<<std::endl
- <<"Time step data:"<<std::endl
- <<" initial_dt = "<<initial_dt<<std::endl
- <<" final_dt = "<<final_dt<<std::endl
- <<" dt_decrement = "<<dt_decrement<<std::endl<<std::endl
- <<"Space discretization data:"<<std::endl
- <<" n_of_global_refines = "<<n_of_global_refines<<std::endl
- <<" Pressure space degree = "<<pressure_degree<<std::endl<<std::endl
- <<"Data for the solution of the diffusion:"<<std::endl
- <<" max iterations = "<<vel_max_iterations<<std::endl
- <<" eps = "<<vel_eps<<std::endl
- <<" Krylov subspace size = "<<vel_Krylov_size<<std::endl
- <<" off diagonals = "<<vel_off_diagonals<<std::endl
- <<" diagonal strengthening = "<<vel_diag_strength<<std::endl
- <<" update prec every "<<vel_update_prec<<" steps"<<std::endl<<std::endl
- <<"Data for the projection step:"<<std::endl
- <<" max iterations = "<<proj_max_iterations<<std::endl
- <<" eps = "<<proj_eps<<std::endl
- <<" off diagonals = "<<proj_off_diagonals<<std::endl
- <<" diagonal strengthening = "<<proj_diag_strength<<std::endl<<std::endl
- <<"Data for the pressure update step"<<std::endl
- <<" max iterations = "<<pres_max_iterations<<std::endl
- <<" eps = "<<pres_eps<<std::endl
- <<" off_diagonals = "<<pres_off_diagonals<<std::endl
- <<" diagonal strengthening = "<<pres_diag_strength<<std::endl<<std::endl
- <<"Verbose = "<<(verbose?"true":"false")<<std::endl<<std::endl
- <<"Output ="<<output<<std::endl<<std::endl;
-}
-
-
-
-// This method just prints the format the
-// input data file should have
-void Data_Storage::print_usage(){
- static const char *message = "";
- std::cout<<message<<std::endl;
- prm.print_parameters( std::cout, ParameterHandler::Text );
-}
+++ /dev/null
-/*
- Muti_Component_Function class implementation
-*/
-#include "../include/MultiComponentFunction.h"
-
-
-// Constructor: Just call the Function constructor and set the component to 0
-template<int dim> Multi_Component_Function<dim>::Multi_Component_Function( const double initial_time ):
- Function<dim>( 1, initial_time ), component(0) {
-}
-
-
-
-// Set Component Function: Check that it is in range and then set it
-template<int dim> void Multi_Component_Function<dim>::set_component(const unsigned int d ){
- // Check if the requested component is correct
- Assert( d<dim, ExcIndexRange( d, 0, dim ) );
- // We have only written 2d
- Assert( dim >= 2, ExcNotImplemented() );
- component = d;
-}
-
-
-
-// explicit template instantiation
-template class Multi_Component_Function<deal_II_dimension>;
+++ /dev/null
-/*
- Implementation of the Navier_Stokes_projection class
- by Abner Salgado.
-*/
-#include "../include/NavierStokes.h"
-
-
-
-#include <sys/time.h>
-
-/// debug
-bool show;
-double Solve_time;
-///----
-
-
-// Constructor
-template<int dim> Navier_Stokes_Projection<dim>::Navier_Stokes_Projection( const 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 ), 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 )
-{
- // After having initialized the bunch of data we do nothing
- // NOTE TO SELF: Do I need to do this check?
- if(deg < 1)
- std::cout<<" WARNING: The chosen pair of finite element spaces is not stable."<<std::endl
- <<" The obtained results will be nonsense"<<std::endl;
- // We check that the time step is within the allowed limits
- AssertThrow( not ( ( dt <= 0. ) or ( dt > .5*T ) ), ExcInvalidTimeStep( dt, .5*T ) );
-}
-
-
-
-// Destructor
-template<int dim> Navier_Stokes_Projection<dim>::~Navier_Stokes_Projection(){
- dof_handler_velocity.clear();
- dof_handler_pressure.clear();
-}
-
-
-
-// Set time step
-template<int dim> void Navier_Stokes_Projection<dim>::set_dt( const double ddt ){
- // We just check that it is within the permitted limits
- AssertThrow( not ( ( ddt <= 0. ) or ( ddt > .5*T ) ), ExcInvalidTimeStep( ddt, .5*T ) );
- dt = ddt;
-}
-
-
-
-// Initialization of the velocity matrices and assembly of those that do not depend on dt
-template<int dim> void Navier_Stokes_Projection<dim>::init_velocity_matrices(){
- //// Init the sparsity pattern for the velocity
- spar_pattern_velocity.reinit( dof_handler_velocity.n_dofs(), dof_handler_velocity.n_dofs(),
- dof_handler_velocity.max_couplings_between_dofs() );
- DoFTools::make_sparsity_pattern( dof_handler_velocity, spar_pattern_velocity );
- spar_pattern_velocity.compress();
-
- //// Init the matrices for the velocity
- vel_Laplace_plus_Mass.reinit( spar_pattern_velocity );
- for( unsigned int d=0; d<dim; ++d )
- vel_it_matrix[d].reinit( spar_pattern_velocity );
- vel_Mass.reinit( spar_pattern_velocity );
- vel_Laplace.reinit( spar_pattern_velocity );
-
- // We assemble the mass matrix
- MatrixCreator::create_mass_matrix( dof_handler_velocity, quadrature_velocity, vel_Mass );
-
- // We assemble the Laplace matrix
- MatrixCreator::create_laplace_matrix( dof_handler_velocity, quadrature_velocity, vel_Laplace );
-}
-
-
-
-template<int dim> void Navier_Stokes_Projection<dim>::init_pressure_matrices(){
- //// Init the sparsity pattern for the pressure
- spar_pattern_pressure.reinit( dof_handler_pressure.n_dofs(), dof_handler_pressure.n_dofs(),
- dof_handler_pressure.max_couplings_between_dofs() );
- DoFTools::make_sparsity_pattern( dof_handler_pressure, spar_pattern_pressure );
-
- // Before we close the sparsity pattern, we need to
- // init the constraints for the Laplace operator on the pressure space
- pres_regularization.clear();
-// DoFTools::make_hanging_node_constraints( dof_handler_pressure, pres_regularization ); //??
-
- // Add the only constraint that we have
- pres_regularization.add_line(0);
-
- // close it
- pres_regularization.close();
-
- // condense the sparsity pattern
- pres_regularization.condense( spar_pattern_pressure );
-
- // Compress the sparsity pattern for the pressure
- spar_pattern_pressure.compress();
-
- //// Init the matrices for the pressure
- pres_Laplace.reinit( spar_pattern_pressure );
- pres_Mass.reinit( spar_pattern_pressure );
-
- // Now we assemble the matrices
- // The Laplace operator is the projection matrix
- MatrixCreator::create_laplace_matrix( dof_handler_pressure, quadrature_pressure, pres_Laplace );
- // The pressure mass matrix to do the pressure update
- MatrixCreator::create_mass_matrix( dof_handler_pressure, quadrature_pressure, pres_Mass );
-
- // Finally we condense the Laplace operator
- pres_regularization.condense( pres_Laplace );
-}
-
-
-
-template<int dim> void Navier_Stokes_Projection<dim>::init_gradient_operator(){
- //// Init the sparsity pattern for the gradient operator
- spar_pattern_pres_vel.reinit( dof_handler_velocity.n_dofs(), dof_handler_pressure.n_dofs(),
- dof_handler_velocity.max_couplings_between_dofs() );
- DoFTools::make_sparsity_pattern( dof_handler_velocity, dof_handler_pressure, spar_pattern_pres_vel );
- spar_pattern_pres_vel.compress();
-
- /*
- To assemble each component of the gradient operator
- we need to make a loop over all cells and compute the products of
- velocity * d_#(pressure)
- where # is the current space component. For this reason we need two cell
- iterators, one for the velocity and one for the pressure space.
- */
- typename DoFHandler<dim>::active_cell_iterator cell_init = dof_handler_velocity.begin_active(),
- cell_end = dof_handler_velocity.end(),
- cell,
- bogus_cell_init = dof_handler_pressure.begin_active(),
- bogus_cell;
- /*
- The FEValues extractors.
- For the velocity we need values,
- For the pressure we only need gradients.
- */
- FEValues<dim> fe_values_velocity( fe_velocity, quadrature_velocity, update_values | update_JxW_values ),
- fe_values_pressure( fe_pressure, quadrature_velocity, update_gradients );
-
- // Usual, and useful, abreviations
- const unsigned int vel_dofs_per_cell = fe_velocity.dofs_per_cell,
- pres_dofs_per_cell = fe_pressure.dofs_per_cell,
- n_q_points = quadrature_velocity.size();
-
- // The local gradient operator
- FullMatrix<double> local_grad( vel_dofs_per_cell, pres_dofs_per_cell );
-
- // Local to global DoF's map
- std::vector<unsigned int> vel_local_dof_indices( vel_dofs_per_cell ), pres_local_dof_indices( pres_dofs_per_cell );
-
- for( unsigned int d=0; d<dim; ++d ){
- // we reinit the given component of the gradient operator
- pres_Diff[d].reinit( spar_pattern_pres_vel );
-
- // The usual cell loop
- for( cell = cell_init, bogus_cell = bogus_cell_init; cell not_eq cell_end; ++cell, ++bogus_cell ){
- fe_values_velocity.reinit( cell );
- fe_values_pressure.reinit( bogus_cell );
-
- cell->get_dof_indices( vel_local_dof_indices );
- bogus_cell->get_dof_indices( pres_local_dof_indices );
-
- // local contributions
- local_grad = 0.;
- for( unsigned int q=0; q<n_q_points; ++q ){
- for( unsigned int i=0; i<vel_dofs_per_cell; ++i )
- for( unsigned int j=0; j<pres_dofs_per_cell; ++j )
- local_grad( i, j ) += fe_values_velocity.JxW(q)*fe_values_velocity.shape_value( i, q )
- *fe_values_pressure.shape_grad( j, q)[d];
- }
-
- // Add it to the global contributions
- for( unsigned int i=0; i<vel_dofs_per_cell; ++i )
- for( unsigned int j=0; j<pres_dofs_per_cell; ++j)
- pres_Diff[d].add( vel_local_dof_indices[i], pres_local_dof_indices[j], local_grad( i, j) );
- }
- }
-}
-
-
-
-// Create Triangulation method
-template<int dim> void Navier_Stokes_Projection<dim>::Create_Triangulation( const unsigned int n_of_refines ){
- // A disk
- GridGenerator::hyper_ball( triangulation );
- static const HyperBallBoundary<dim> boundary;
- triangulation.set_boundary( 0, boundary );
-
-/* // Our domain is a unit square.
- GridGenerator::hyper_cube(triangulation);*/
- triangulation.refine_global( n_of_refines );
- std::cout<<" Number of active cells: "<<triangulation.n_active_cells()<<std::endl;
-
- // Having created the mesh, we create the mesh-dependent data
-
- //// Distribute the degrees of freedom and renumber them
- dof_handler_velocity.distribute_dofs( fe_velocity );
- DoFRenumbering::boost::Cuthill_McKee( dof_handler_velocity );
- dof_handler_pressure.distribute_dofs( fe_pressure );
- DoFRenumbering::boost::Cuthill_McKee( dof_handler_pressure );
-
- // Initialize the matrices for the velocity
- init_velocity_matrices();
-
- // Initialize the matrices for the pressure
- init_pressure_matrices();
-
- // Initialize the gradient operator
- init_gradient_operator();
-
-
- //// Set the correct size for the various vectors that we need
- 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() );
- phi_n_minus_1.reinit( dof_handler_pressure.n_dofs() );
-
- pres_tmp.reinit( dof_handler_pressure.n_dofs() );
-
- for(unsigned int d=0; d<dim; ++d){
- u_n[d].reinit( dof_handler_velocity.n_dofs() );
- u_n_minus_1[d].reinit( dof_handler_velocity.n_dofs() );
- u_star[d].reinit( dof_handler_velocity.n_dofs() );
- force[d].reinit( dof_handler_velocity.n_dofs() );
- }
- v_tmp.reinit( dof_handler_velocity.n_dofs() );
-
- std::cout<<" dim( X_h ) = "<<( dof_handler_velocity.n_dofs()*dim )<<std::endl
- <<" dim( M_h ) = "<<dof_handler_pressure.n_dofs()<<std::endl;
-}
-
-
-
-// Initialization of the time dependent data
-template<int dim> void Navier_Stokes_Projection<dim>::Initialize(){
-
- // The constant part of the matrix for the velocity is 1.5/dt*Mass + 1/Re*Delta_h
- vel_Laplace_plus_Mass = 0.;
- vel_Laplace_plus_Mass.add( 1./Re, vel_Laplace );
- vel_Laplace_plus_Mass.add( 1.5/dt, vel_Mass );
-
- // Now we start the vectors and load them with the initial data
- ////pressure
- Pressure<dim> pres( t_0 );
- VectorTools::interpolate( dof_handler_pressure, pres, pres_n_minus_1 );
- pres.advance_time( dt );
- VectorTools::interpolate( dof_handler_pressure, pres, pres_n );
-
- ////phi
- phi_n = 0.;
- phi_n_minus_1 = 0.;
-
- ////velocity
- Velocity<dim> v;
- for(unsigned int d=0; d<dim; ++d){
- v.set_time( t_0 );
- v.set_component(d);
-
- VectorTools::interpolate( dof_handler_velocity, v, u_n_minus_1[d] );
-
- v.advance_time( dt );
- VectorTools::interpolate( dof_handler_velocity, v, u_n[d] );
- }
- plot_solution(1);
-}
-
-
-
-/*
- A time marching step, the output messages and name of the functions are self
- explanatory, so no comments here
-*/
-template<int dim> void Navier_Stokes_Projection<dim>::run( const bool verbose, const unsigned int n_of_plots ){
- // We set the number of steps
- unsigned int n_steps = ( T - t_0 )/dt;
- // Set the source term to the correct time
- rhs.set_time( 2.*dt );
- vel_exact.set_time( 2.*dt );
- // do as many steps as necessary
- timeval init_time, end_time;
-
- /// debug
- show = false;
- Solve_time = 0.;
- ///-----
-
- gettimeofday( &init_time, 0 );
- for( unsigned int n = 2; n<=n_steps; ++n ){
-
- /// debug
- if( n == n_steps )
- show = true;
- ///----
-
- if( verbose )
- std::cout<<" Step = "<<n<<" Time = "<<(n*dt)<<std::endl;
- if( verbose )
- std::cout<<" Interpolating the velocity "<<std::endl;
- interpolate_velocity();
- if( verbose )
- std::cout<<" Diffusion Step"<<std::endl;
- if( ( n%vel_update_prec==0 ) and ( verbose ) )
- std::cout<<" With reinitialization of the preconditioner"<<std::endl;
- diffusion_step( (n%vel_update_prec == 0 ) or ( n == 2) );
- if( verbose )
- std::cout<<" Projection Step"<<std::endl;
- projection_step( ( n == 2 ) );
- if( verbose )
- std::cout<<" Updating the Pressure"<<std::endl;
- update_pressure( ( n == 2 ) );
- if( n%n_of_plots ==0 ){
- if( verbose )
- std::cout<<" Plotting Solution"<<std::endl;
- plot_solution(n);
- }
- // advance in time
- rhs.advance_time(dt);
- vel_exact.advance_time(dt);
- }
- gettimeofday( &end_time, 0 );
- double run_time = double( end_time.tv_sec - init_time.tv_sec ) + 1e-6*double( end_time.tv_usec - init_time.tv_usec );
- run_time /= n_steps;
- run_time /= double( dim*dof_handler_velocity.n_dofs() + dof_handler_pressure.n_dofs() );
-// if( verbose )
- std::cout<<" Run time = "<<run_time<<std::endl;
-}
-
-
-
-// Simple procedure that interpolates the velocity
-template<int dim> void Navier_Stokes_Projection<dim>::interpolate_velocity(){
- for( unsigned int d=0; d<dim; ++d ){
- u_star[d] = 0.;
- u_star[d].add( 2., u_n[d], -1., u_n_minus_1[d] );
- }
-}
-
-
-
-template<int dim> void Navier_Stokes_Projection<dim>::assemble_advection_term( const unsigned int d ){
- // Cell iterator
- typename DoFHandler<dim>::active_cell_iterator cell = dof_handler_velocity.begin_active(),
- cell_end = dof_handler_velocity.end();
-
- /*
- The FEValues extractors.
- We need values because of the terms involving u^*
- also we need the gradient to assemble the advection part
- */
- FEValues<dim> fe_values( fe_velocity, quadrature_velocity, update_values | update_JxW_values | update_gradients );
-
- // Usual, and useful, abreviations
- const unsigned int dofs_per_cell = fe_velocity.dofs_per_cell,
- n_q_points = quadrature_velocity.size();
-
- /*
- We need to compute the values of several finite element functions (and some gradients)
- at quadrature points. For this we use all these arrays.
- The names are self-explanatory
- */
- std::vector<double> u_old_local( n_q_points );
- std::vector< Point<dim> > u_star_local( n_q_points );
- std::vector< Tensor<1,dim> > grad_u_star(n_q_points);
-
- // Local to global DoF's map
- std::vector<unsigned int> local_dof_indices( dofs_per_cell );
-
- // Local contribution of the advection operator
- Vector<double> local_rhs( dofs_per_cell );
-
- // loop over cells
- for( ; cell not_eq cell_end; ++cell ){
- // reinit the stuff
- fe_values.reinit( cell );
- cell->get_dof_indices( local_dof_indices );
-
- // create the local u^*
- for( unsigned int s=0; s<dim; ++s ){
- fe_values.get_function_values( u_star[s], u_old_local );
- for( unsigned int q=0; q<n_q_points; ++q )
- u_star_local[q](s) = u_old_local[q];
- }
-
- fe_values.get_function_gradients( u_star[d], grad_u_star );
-
- local_rhs = 0.;
- for( unsigned int q =0; q<n_q_points; ++q )
- for( unsigned int i=0; i<dofs_per_cell; ++i )
- // The local contribution of the advection term
- local_rhs(i) -= ( ( u_star_local[q]*grad_u_star[q] )*fe_values.shape_value( i, q ) )*fe_values.JxW(q) ;
-
- // Add local contributions to the global matrix
- for( unsigned int i=0; i<dofs_per_cell; ++i )
- force[d]( local_dof_indices[i] ) += local_rhs(i);
- }
-}
-
-
-
-
-template<int dim> void Navier_Stokes_Projection<dim>::threaded_solve( const unsigned int d){
-
- // Create the solver that we are going to use
- SolverControl solver_control( vel_max_its, vel_eps*force[d].l2_norm() );
- SolverCG<> cg( solver_control );
- // Solve
- cg.solve( vel_it_matrix[d], u_n[d], force[d], prec_velocity[d] );
-}
-
-
-
-// A diffusion step
-template<int dim> void Navier_Stokes_Projection<dim>::diffusion_step( const bool reinit_prec ){
- // Used to solve in parallel the components of the velocity
- Threads::ThreadGroup<> threads;
-
- /// debug
- static double assembly_time = 0.;
- timeval init_time_assembly, end_time_assembly, init_time_solve, end_time_solve;
- gettimeofday( &init_time_assembly, 0 );
- ///----
-
- /*
- Define the pressure interpolant. It is defined by
- p^* = p^{n} + 4/3 phi^{n} - 1/3 phi^{n-1}
- */
- pres_tmp = pres_n;
- pres_tmp.add(4./3., phi_n, -1./3., phi_n_minus_1);
- pres_tmp *= -1.;
-
- for( unsigned int d=0; d<dim; ++d ){
- // First we need to construct the right hand sides
- rhs.set_component(d);
- //// compute the source term
- VectorTools::create_right_hand_side( dof_handler_velocity, quadrature_velocity, rhs, force[d] );
-
- // We assemble the advection term
- threads += Threads::spawn( *this, &Navier_Stokes_Projection<dim>::assemble_advection_term )(d);
- //assemble_advection_term(d);
-
- /*
- Compute the value of the contribution of the old velocities.
- Recall that BDF2 reads
- ( 0.5/dt)*( 3u^{n+1} - 4u^{n} + u^{n-1} )
- the terms involving u^{n} and u^{n-1} are known, so they go into the
- right hand side
- */
- v_tmp = 0.;
- v_tmp.add( 2./dt,u_n[d],-.5/dt,u_n_minus_1[d] );
- vel_Mass.vmult_add( force[d], v_tmp );
-
- // And add the contribution of the pressure part
- pres_Diff[d].vmult_add( force[d], pres_tmp );
-
- // Copy the old value
- u_n_minus_1[d] = u_n[d];
-
-
- // Create the matrix
- vel_it_matrix[d].copy_from( vel_Laplace_plus_Mass );
- //// Create the boundary values
- vel_exact.set_component(d);
- VectorTools::interpolate_boundary_values( dof_handler_velocity, 0, vel_exact, boundary_values );
-
- // apply BCs
- MatrixTools::apply_boundary_values( boundary_values, vel_it_matrix[d], u_n[d], force[d] );
- }
- threads.join_all();
-
- /// debug
- gettimeofday( &end_time_assembly, 0 );
- assembly_time += double( end_time_assembly.tv_sec - init_time_assembly.tv_sec )
- + 1e-6*double( end_time_assembly.tv_usec - init_time_assembly.tv_usec );
- gettimeofday( &init_time_solve, 0 );
- ///-----
-
- for(unsigned int d=0; d<dim; ++d ){
- // Check if the preconditioner needs reinitialization
- // if( reinit_prec )
- prec_velocity[d].initialize( vel_it_matrix[d]/*, SparseILU<double>::AdditionalData( vel_diag_strength, vel_off_diagonals )*/ );
-
- //prec_velocity[d].factorize( vel_it_matrix[d] );
-// threaded_solve(d);
-
-// And solve
-// threads += Threads::spawn( *this, &Navier_Stokes_Projection<dim>::threaded_solve )(d);
-
-// prec_velocity[d].solve( force[d] );
-// u_n[d] = force[d];
- prec_velocity[d].vmult( force[d], u_n[d] );
- }
-// threads.join_all();
-
-
- ///debug
- gettimeofday( &end_time_solve, 0 );
- Solve_time += double( end_time_solve.tv_sec - init_time_solve.tv_sec )
- + 1e-6*double( end_time_solve.tv_usec - init_time_solve.tv_usec );
- if(show)
- std::cout<<" Assembly time = "<<assembly_time<<std::endl
- <<" Solve time = "<<Solve_time<<std::endl;
- ///----
-}
-
-
-
-// A projection step
-template<int dim> void Navier_Stokes_Projection<dim>::projection_step( const bool reinit_prec ){
-
- /// debug
- timeval init_time_assemble, init_time_solve, end_time_assemble, end_time_solve, init_time_factorize, end_time_factorize;
- static double assemble_time=0., solve_time=0., factorize_time = 0.;
- ///----
-
-
- // Check if the preconditioner needs reinitialization
- if( reinit_prec ){
-
- /// debug
- gettimeofday( &init_time_factorize, 0 );
- ///-----
-// prec_pressure.initialize( pres_Laplace, SparseILU<double>::AdditionalData( proj_diag_strength, proj_off_diagonals ) );
- prec_pressure.initialize( pres_Laplace );
-
- /// debug
- gettimeofday( &end_time_factorize, 0 );
- factorize_time += double( end_time_factorize.tv_sec - init_time_factorize.tv_sec )
- + 1e-6*double( end_time_factorize.tv_usec - init_time_factorize.tv_usec );
- ///--------
- }
-
- /// debug
- gettimeofday( &init_time_assemble, 0 );
- ///-----
-
- // To make things faster we assembled the discrete gradient operator
- // so we just need to multiply by its transpose
- pres_tmp = 0.;
- for( unsigned d=0; d<dim; ++d )
- pres_Diff[d].Tvmult_add( pres_tmp, u_n[d] );
-
- /// debug
- gettimeofday( &end_time_assemble, 0 );
- assemble_time += double( end_time_assemble.tv_sec - init_time_assemble.tv_sec )
- + 1e-6*double( end_time_assemble.tv_usec - init_time_assemble.tv_usec );
- ///-----
-
- // Up to this point, pres_tmp has the term < u_n, grad q >
- // we will assemble this only once, so we
- // Copy the old value of phi
- phi_n_minus_1 = phi_n;
-
- phi_n = pres_tmp;
- phi_n *= 1.5/dt;
-
- /// debug
- gettimeofday( &init_time_solve, 0 );
- ///-------
-
- // condense because of the constraint
- pres_regularization.condense( phi_n );
-
-
- prec_pressure.solve( phi_n );
- // The method that we are going to use
-// SolverControl solver_control( proj_max_its, proj_eps*pres_tmp.l2_norm() );
-// SolverCG<> cg( solver_control );
-
-
- // And solve
-// cg.solve( pres_Laplace, phi_n, pres_tmp, prec_pressure );
-
- // After solving we need to distribute the constrained stuff
- pres_regularization.distribute( phi_n );
-
- /// debug
- gettimeofday( &end_time_solve, 0 );
- solve_time += double( end_time_solve.tv_sec - init_time_solve.tv_sec )
- + 1e-6*double( end_time_solve.tv_usec - init_time_solve.tv_usec );
- if(show)
- std::cout<<" DIV assembly time = "<<assemble_time<<std::endl
- <<" LAP solve time = "<<solve_time<<std::endl
- <<" LAP factorize time = "<<factorize_time<<std::endl;
- ///------
-}
-
-
-
-// Pressure update step
-template<int dim> void Navier_Stokes_Projection<dim>::update_pressure( const bool reinit_prec ){
-
- /// debug
- static double solve_time = 0.;
- timeval init_time_solve, end_time_solve;
- ///---------
-
- // Copy the old pressure
- pres_n_minus_1 = pres_n;
- switch( type ){
- case METHOD_STANDARD:
- /*
- Simple standard case
- p^{n+1} = p^n + \phi^{n+1}
- */
- pres_n += phi_n;
- break;
- case METHOD_ROTATIONAL:
- {
- /*
- Rotational form
- p^{n+1} = p^n + \phi^{n+1} - 1/Re div u^{n+1}
- */
-
- /// debug
- gettimeofday( &init_time_solve, 0 );
- ///------
-
- // The preconditioner for the mass matrix problem
- if( reinit_prec )
- prec_mass.initialize( pres_Mass/*, SparseILU<double>::AdditionalData( pres_diag_strength, pres_off_diagonals )*/ );
-
- // The solver for the mass matrix
-// SolverControl solver_control( pres_max_its, pres_eps*pres_tmp.l2_norm() );
-// SolverCG<> cg( solver_control );
-//
-// // Solve the mass matrix problem
-// cg.solve( pres_Mass, pres_n, pres_tmp, prec_mass );
-
- pres_n = pres_tmp;
- prec_mass.solve( pres_n );
-
- pres_n.sadd(1./Re, 1., pres_n_minus_1, 1., phi_n );
-
- /// debug
- gettimeofday( &end_time_solve, 0 );
- solve_time += double( end_time_solve.tv_sec - init_time_solve.tv_sec )
- + 1e-6*double( end_time_solve.tv_usec - init_time_solve.tv_usec );
- ///----------
-
- break;
- }
- default:
- Assert( false, ExcNotImplemented() );
- };
-
- /// debug
- if(show)
- std::cout<<" PRES_UPD solve time = "<<solve_time<<std::endl;
- ///---
-
-}
-
-
-
-// Post Processing
-template<int dim> void Navier_Stokes_Projection<dim>::Post_Process(){
- double tmp, vel_err_L2=0., vel_err_H1=0., pres_err_L2;
-
- Vector<double> differences( triangulation.n_active_cells() );
-
- Velocity<dim> u_exact(T);
- for( unsigned int d=0; d<dim; ++d ){
- u_exact.set_component(d);
-
- // L2 error for this component of the velocity
- differences = 0.;
- VectorTools::integrate_difference( dof_handler_velocity, u_n[d], u_exact, differences, quadrature_velocity, VectorTools::L2_norm );
- tmp = differences.l2_norm();
- vel_err_L2 += tmp*tmp;
-
- // H1 error for this component of the velocity
- differences = 0.;
- VectorTools::integrate_difference( dof_handler_velocity, u_n[d], u_exact, differences, quadrature_velocity, VectorTools::H1_seminorm );
- tmp = differences.l2_norm();
- vel_err_H1 += tmp*tmp;
- }
- vel_err_L2 = std::sqrt( vel_err_L2 );
- vel_err_H1 = std::sqrt( vel_err_H1 );
-
-
- /*
- The actual pressure needs to have mean value zero.
- We compute the mean value and subtract it from the computed pressure
- */
- double pres_mean_value = VectorTools::compute_mean_value( dof_handler_pressure, quadrature_pressure, pres_n, 0 );
- for( unsigned int i=0; i<pres_n.size(); ++i )
- pres_n(i) -= pres_mean_value;
-
- Pressure<dim> pres_exact(T);
-
- // L2 error of the pressure
- differences = 0.;
- VectorTools::integrate_difference( dof_handler_pressure, pres_n, pres_exact, differences, quadrature_pressure, VectorTools::L2_norm );
- pres_err_L2 = differences.l2_norm();
-
- // Add the obtained results to the convergence table
- 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 );
-
- // And show it
- convergence_table.write_text(std::cout);
-}
-
-
-
-// Output of the solution
-template<int dim> void Navier_Stokes_Projection<dim>::plot_solution( const unsigned int step ){
- // This only works in 2d
- Assert( dim==2, ExcNotImplemented() );
-
- // We need to assemble the vorticy
- FEValues<dim> fe_values( fe_velocity, quadrature_velocity, update_values | update_gradients | update_JxW_values );
- const unsigned int dofs_per_cell = fe_velocity.n_dofs_per_cell(),
- n_q_points = quadrature_velocity.size();
- std::vector< Tensor<1,dim> > grad_vel_1( n_q_points ), grad_vel_2( n_q_points );
- std::vector< unsigned int> local_dof_indices( dofs_per_cell );
- Vector<double> local_rhs( dofs_per_cell );
- double vorticity;
-
- typename DoFHandler<dim>::active_cell_iterator cell = dof_handler_velocity.begin_active(),
- cend = dof_handler_velocity.end();
-
- force[0] = 0.;
- //We start the usual loop
- for( ; cell not_eq cend; ++cell ){
- // reinit all the needed stuff
- local_rhs = 0.;
- fe_values.reinit( cell );
- cell->get_dof_indices( local_dof_indices );
-
- // get the gradients of each function
- fe_values.get_function_gradients( u_n[0], grad_vel_1 );
- fe_values.get_function_gradients( u_n[1], grad_vel_2 );
-
- // usual loop over quad points and local dofs
- for( unsigned int q=0; q<n_q_points; ++q ){
- vorticity = grad_vel_2[q][1] - grad_vel_1[q][0];
- for( unsigned int i=0; i<dofs_per_cell; ++i )
- local_rhs(i) += fe_values.JxW(q)*vorticity*fe_values.shape_value( i, q );
- }
-
- // And copy it to the local one
- for( unsigned int i=0; i<dofs_per_cell; ++i )
- force[0]( local_dof_indices[i] ) += local_rhs(i);
- }
-
- // We solve a mass matrix problem
- SolverControl solver_control( 500, 1e-6*force[0].l2_norm() );
- SolverCG<> cg( solver_control );
- static bool is_prec_initted = false;
- static SparseILU<double> prec;
- if( not is_prec_initted ){
- prec.initialize( vel_Mass, SparseILU<double>::AdditionalData( 1e-5, 70 ) );
- is_prec_initted = true;
- }
- cg.solve( vel_Mass, force[1], force[0], prec );
-
- // Once we have the vorticity we can output it
- DataOut<dim> data_out;
- data_out.attach_dof_handler( dof_handler_velocity );
- data_out.add_data_vector( force[1], "vorticity" );
- data_out.build_patches();
- std::ostringstream filename;
- filename<<"vorticity"<<step<<".vtk";
- std::ofstream output( filename.str().c_str() );
- data_out.write_vtk( output );
-}
-
-
-// explicit template instantiation
-template class Navier_Stokes_Projection<deal_II_dimension>;
+++ /dev/null
-#include "../include/NavierStokes.h"
-
-
-#include <iostream>
-#include <sys/time.h>
-
-
-
-int main( int argc, char **argv ){
- try{
- Data_Storage data;
- if( argc<2 ){
- std::cout<<std::endl<<std::endl
- <<" Missing parameter file!!"<<std::endl
- <<" Aborting."<<std::endl<<std::endl;
- data.print_usage();
- exit(1);
- }
- data.read_data( argv[1] );
- deallog.depth_console( data.verbose?2:0 );
- Navier_Stokes_Projection<deal_II_dimension> test( data );
- test.Create_Triangulation( data.n_of_global_refines );
- timeval init_time, end_time;
- gettimeofday( &init_time, 0 );
- for( double dt = data.initial_dt; dt >= data.final_dt; dt /= data.dt_decrement ){
- std::cout<<" dt = "<<dt<<std::endl;
- test.set_dt( dt );
- test.Initialize();
- test.run( data.verbose, data.output );
- test.Post_Process();
- std::cout<<"====================================="<<std::endl<<std::endl;
- }
- gettimeofday( &end_time, 0 );
- std::cout<<" Global time = "
- <<( double( end_time.tv_sec - init_time.tv_sec )+1e-6*double( end_time.tv_usec - init_time.tv_usec ) )
- <<" sec"<<std::endl;
- }
- catch (std::exception &exc){
- std::cerr << std::endl << std::endl
- << "----------------------------------------------------"
- << std::endl;
- std::cerr << "Exception on processing: " << std::endl
- << exc.what() << std::endl
- << "Aborting!" << std::endl
- << "----------------------------------------------------"
- << std::endl;
- return 1;
- }
- catch (...){
- std::cerr << std::endl << std::endl
- << "----------------------------------------------------"
- << std::endl;
- std::cerr << "Unknown exception!" << std::endl
- << "Aborting!" << std::endl
- << "----------------------------------------------------"
- << std::endl;
- return 1;
- }
- std::cout<<"----------------------------------------------------"
- <<std::endl
- <<"Apparently everything went fine!"
- <<std::endl
- <<"Don't forget to brush your teeth :-)"
- <<std::endl<<std::endl;
- return 0;
-}
projects / dealii-svn.git / commitdiff