#include <fe/fe_q.h>
#include <fe/fe_dgq.h>
+#include <fe/fe_dgp.h>
#include <fe/fe_system.h>
#include <fe/fe_values.h>
#include <fe/mapping_q.h>
unsigned int adaptive_refinement_interval;
- double stabilization_alpha;
- double stabilization_c_R;
- double stabilization_beta;
+ double stabilization_alpha;
+ double stabilization_c_R;
+ double stabilization_beta;
+
+ unsigned int stokes_velocity_degree;
+ bool use_locally_conservative_discretization;
+
+ unsigned int temperature_degree;
};
const MappingQ<dim> mapping;
- const unsigned int stokes_degree;
const FESystem<dim> stokes_fe;
DoFHandler<dim> stokes_dof_handler;
TrilinosWrappers::MPI::BlockVector stokes_rhs;
- const unsigned int temperature_degree;
FE_Q<dim> temperature_fe;
DoFHandler<dim> temperature_dof_handler;
ConstraintMatrix temperature_constraints;
adaptive_refinement_interval (10),
stabilization_alpha (2),
stabilization_c_R (0.11),
- stabilization_beta (0.078)
+ stabilization_beta (0.078),
+ stokes_velocity_degree (2),
+ use_locally_conservative_discretization (true),
+ temperature_degree (2)
{}
"and 0.078 for 3d.");
}
prm.leave_subsection ();
+
+ prm.enter_subsection ("Discretization");
+ {
+ prm.declare_entry ("Stokes velocity polynomial degree", "2",
+ Patterns::Integer (1),
+ "The polynomial degree to use for the velocity variables "
+ "in the Stokes system.");
+ prm.declare_entry ("Temperature polynomial degree", "2",
+ Patterns::Integer (1),
+ "The polynomial degree to use for the temperature variable.");
+ prm.declare_entry ("Use locally conservative discretization", "true",
+ Patterns::Bool (),
+ "Whether to use a Stokes discretization that is locally "
+ "conservative at the expense of a larger number of degrees "
+ "of freedom, or to go with a cheaper discretization "
+ "that does not locally conserve mass (although it is "
+ "globally conservative.");
+ }
+ prm.leave_subsection ();
}
stabilization_beta = prm.get_double ("beta");
}
prm.leave_subsection ();
+
+ prm.enter_subsection ("Discretization");
+ {
+ stokes_velocity_degree = prm.get_integer ("Stokes velocity polynomial degree");
+ temperature_degree = prm.get_integer ("Temperature polynomial degree");
+ use_locally_conservative_discretization
+ = prm.get_bool ("Use locally conservative discretization");
+ }
+ prm.leave_subsection ();
}
mapping (4),
- stokes_degree (1),
- stokes_fe (FE_Q<dim>(stokes_degree+1), dim,
- FE_Q<dim>(stokes_degree), 1),
+ stokes_fe (FE_Q<dim>(parameters.stokes_velocity_degree),
+ dim,
+ (parameters.use_locally_conservative_discretization
+ ?
+ static_cast<const FiniteElement<dim> &>
+ (FE_DGP<dim>(parameters.stokes_velocity_degree-1))
+ :
+ static_cast<const FiniteElement<dim> &>
+ (FE_Q<dim>(parameters.stokes_velocity_degree-1))),
+ 1),
stokes_dof_handler (triangulation),
- temperature_degree (2),
- temperature_fe (temperature_degree),
+ temperature_fe (parameters.temperature_degree),
temperature_dof_handler (triangulation),
time_step (0),
double BoussinesqFlowProblem<dim>::get_maximal_velocity () const
{
const QIterated<dim> quadrature_formula (QTrapez<1>(),
- stokes_degree+1);
+ parameters.stokes_velocity_degree);
const unsigned int n_q_points = quadrature_formula.size();
FEValues<dim> fe_values (mapping, stokes_fe, quadrature_formula, update_values);
BoussinesqFlowProblem<dim>::get_extrapolated_temperature_range () const
{
const QIterated<dim> quadrature_formula (QTrapez<1>(),
- temperature_degree);
+ parameters.temperature_degree);
const unsigned int n_q_points = quadrature_formula.size();
FEValues<dim> fe_values (mapping, temperature_fe, quadrature_formula,
{
assemble_temperature_matrix ();
- QGauss<dim> quadrature(temperature_degree+2);
+ QGauss<dim> quadrature(parameters.temperature_degree+2);
UpdateFlags update_flags = UpdateFlags(update_values |
update_quadrature_points |
update_JxW_values);
{
stokes_preconditioner_matrix = 0;
- const QGauss<dim> quadrature_formula(stokes_degree+2);
+ const QGauss<dim> quadrature_formula(parameters.stokes_velocity_degree+1);
typedef
FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
stokes_rhs=0;
- const QGauss<dim> quadrature_formula(stokes_degree+2);
+ const QGauss<dim> quadrature_formula(parameters.stokes_velocity_degree+1);
typedef
FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
temperature_mass_matrix = 0;
temperature_stiffness_matrix = 0;
- const QGauss<dim> quadrature_formula(temperature_degree+2);
+ const QGauss<dim> quadrature_formula(parameters.temperature_degree+2);
typedef
FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
temperature_rhs = 0;
- const QGauss<dim> quadrature_formula(temperature_degree+2);
+ const QGauss<dim> quadrature_formula(parameters.temperature_degree+2);
const std::pair<double,double>
global_T_range = get_extrapolated_temperature_range();
// size in 3d
double scaling = (dim==3)?0.25:1.0;
double local_time_step = scaling/(1.6*dim*std::sqrt(1.*dim)) /
- temperature_degree *
+ parameters.temperature_degree *
GridTools::minimal_cell_diameter(triangulation) /
std::max(1e-10,maximal_velocity);
//norm of gradient
{
double my_cells_error[2] = {0, 0};
- QGauss<1> q_base(stokes_degree+1);
+ QGauss<1> q_base(parameters.stokes_velocity_degree);
QIterated<dim> err_quadrature(q_base, 2);
const unsigned int n_q_points = err_quadrature.size();
Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
KellyErrorEstimator<dim>::estimate (temperature_dof_handler,
- QGauss<dim-1>(temperature_degree+1),
+ QGauss<dim-1>(parameters.temperature_degree+1),
typename FunctionMap<dim>::type(),
temperature_solution,
estimated_error_per_cell,
{
parameter_file.close ();
- std::cerr << "*** Input parameter file <"
- << parameter_filename << "> not found. Creating a"
- << std::endl
- << "*** template file of the same name."
- << std::endl;
+ std::ostringstream message;
+ message << "Input parameter file <"
+ << parameter_filename << "> not found. Creating a"
+ << std::endl
+ << "template file of the same name."
+ << std::endl;
std::ofstream parameter_out (parameter_filename.c_str());
prm.print_parameters (parameter_out,
ParameterHandler::Text);
- std::exit (1);
+ AssertThrow (false, ExcMessage (message.str().c_str()));
}
const bool success = prm.read_input (parameter_file);
// This is copied verbatim from step-31:
int main (int argc, char *argv[])
{
+ Utilities::System::MPI_InitFinalize mpi_initialization(argc, argv);
+
try
{
deallog.depth_console (0);
- Utilities::System::MPI_InitFinalize mpi_initialization(argc, argv);
-
std::string parameter_filename;
if (argc>=2)
parameter_filename = argv[1];
# The end time of the simulation in years.
set End time = 1e8
+# Whether graphical output is to be generated or not. You may not want to get
+# graphical output if the number of processors is large.
+set Generate graphical output = false
+
# The number of adaptive refinement steps performed after initial global
# refinement.
set Initial adaptive refinement = 2
# The number of global refinement steps performed on the initial coarse mesh,
# before the problem is first solved there.
-set Initial global refinement = 1
+set Initial global refinement = 2
+
+# The number of time steps between each generation of graphical output files.
+set Time steps between graphical output = 50
# The number of time steps after which the mesh is to be adapted based on
# computed error indicators.
set Time steps between mesh refinement = 10
-# Whether graphical output is to be generated or not. You may not want to get
-# graphical output if the number of processors is large.
-set Generate graphical output = false
-# The number of time steps between each generation of graphical output files.
-set Time steps between graphical output = 50
+subsection Discretization
+ # The polynomial degree to use for the velocity variables in the Stokes
+ # system.
+ set Stokes velocity polynomial degree = 2
+
+ # The polynomial degree to use for the temperature variable.
+ set Temperature polynomial degree = 2
+
+ # Whether to use a Stokes discretization that is locally conservative at the
+ # expense of a larger number of degrees of freedom, or to go with a cheaper
+ # discretization that does not locally conserve mass (although it is
+ # globally conservative.
+ set Use locally conservative discretization = true
+end
subsection Stabilization parameters
# The exponent in the entropy viscosity stabilization.
set alpha = 2
- # The c_R factor in the entropy viscosity stabilization.
- set c_R = 0.11
-
# The beta factor in the artificial viscosity stabilization. An appropriate
# value for 2d is 0.052 and 0.078 for 3d.
set beta = 0.078
+
+ # The c_R factor in the entropy viscosity stabilization.
+ set c_R = 0.11
end