#include <base/timer.h>
-
+#include <base/exceptions.h>
+#include <iostream>
// these includes should probably be properly
// ./configure'd using the AC_HEADER_TIME macro:
#include <sys/time.h>
#include <sys/resource.h>
+#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
+#include <mpi.h>
+#endif
+
// on SunOS 4.x, getrusage is stated in the man pages and exists, but
// is not declared in resource.h. declare it ourselves
running = false;
}
+
+
+/* ---------------------------- TimerOutput -------------------------- */
+
+TimerOutput::TimerOutput (std::ostream &stream,
+ const enum OutputFrequency output_frequency,
+ const enum OutputType output_type)
+ :
+ output_frequency (output_frequency),
+ output_type (output_type),
+ out_stream (stream, true)
+{}
+
+
+
+TimerOutput::TimerOutput (ConditionalOStream &stream,
+ const enum OutputFrequency output_frequency,
+ const enum OutputType output_type)
+ :
+ output_frequency (output_frequency),
+ output_type (output_type),
+ out_stream (stream)
+{}
+
+
+
+TimerOutput::~TimerOutput()
+{
+ while (active_sections.size() > 0)
+ exit_section();
+
+ if (output_frequency != every_call)
+ print_summary();
+}
+
+
+
+void
+TimerOutput::enter_section (const std::string §ion_name)
+{
+ Assert (section_name.empty() == false,
+ ExcMessage ("Section string is empty."));
+
+ unsigned int this_section_number = numbers::invalid_unsigned_int;
+
+ // check whether the requested section
+ // already exists
+ for (unsigned int i=0; i<section_names.size(); ++i)
+ {
+ // if we found it, we have to have the
+ // index into which we write the
+ // data. ensure that the section was not
+ // active before to avoid any mess
+ if (section_names[i].compare (section_name) == 0)
+ {
+ this_section_number = i;
+ for (unsigned int j=0; j<active_sections.size(); ++j)
+ {
+ if (active_sections[j] == i)
+ {
+#ifdef DEBUG
+ std::string exc_text = "Cannot enter the already active section "
+ + section_name;
+ Assert (active_sections[j] != i, ExcMessage (exc_text.data()));
+#endif
+ exit_section(section_name);
+ }
+ }
+ break;
+ }
+ }
+
+ // in case this is a new section, enlarge
+ // our data storage by one, otherwise
+ // just reset the timer
+ if (this_section_number == numbers::invalid_unsigned_int)
+ {
+ section_names.push_back(section_name);
+ section_total_cpu_times.push_back (0.);
+ section_total_wall_times.push_back (0.);
+ section_timers.push_back (Timer());
+ section_n_calls.push_back (0);
+ this_section_number = section_names.size() - 1;
+ }
+ else
+ {
+ section_timers[this_section_number].reset();
+ section_timers[this_section_number].start();
+ }
+
+ // now turn the current section into the
+ // list of active sections and increase
+ // the counter for that section
+ active_sections.push_back (this_section_number);
+ section_n_calls[this_section_number]++;
+}
+
+
+
+void
+TimerOutput::exit_section (const std::string §ion_name)
+{
+ // if no string is given, exit the last
+ // active section.
+ unsigned int active_index_to_delete = numbers::invalid_unsigned_int;
+ if (section_name.empty() == true)
+ active_index_to_delete = active_sections.size() - 1;
+ // if we got a string, we need to find
+ // the index in the list of active
+ // functions
+ else
+ {
+ for (unsigned int i=0; i<active_sections.size(); ++i)
+ if (section_names[active_sections[i]].compare (section_name) == 0)
+ {
+ active_index_to_delete = i;
+ break;
+ }
+ }
+
+ Assert (active_index_to_delete != numbers::invalid_unsigned_int,
+ ExcMessage ("Cannot delete a section that was never created."));
+
+ const unsigned int section_to_exit = active_sections[active_index_to_delete];
+
+ Assert (section_to_exit < section_names.size(),
+ ExcInternalError());
+
+ // stop the timer for this section.
+ section_timers[section_to_exit].stop();
+ section_total_wall_times[section_to_exit] += section_timers[section_to_exit].wall_time();
+
+ // get cpu time. on MPI systems, add the
+ // local contributions.
+ //
+ // TODO: this should rather be in the Timer
+ // class itself, shouldn't it?
+ double cpu_time = section_timers[section_to_exit]();
+ {
+
+ // On MPI, sum up all the local CPU
+ // times.
+#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
+ double total_cpu_time = 0.;
+ int mpiInitialized;
+ MPI_Initialized(&mpiInitialized);
+
+ if( mpiInitialized )
+ {
+ MPI_Allreduce (&cpu_time, &total_cpu_time, 1, MPI_DOUBLE, MPI_SUM,
+ MPI_COMM_WORLD);
+ cpu_time = total_cpu_time;
+ }
+#endif
+ section_total_cpu_times[section_to_exit] += cpu_time;
+ }
+
+
+ // in case we have to print out
+ // something, do that here...
+ if (output_frequency != summary)
+ {
+ std::string output_time;
+ std::ostringstream cpu;
+ cpu << cpu_time << "s";
+ std::ostringstream wall;
+ wall << section_timers[section_to_exit].wall_time() << "s";
+ if (output_type == cpu_times)
+ output_time = " CPU time: " + cpu.str();
+ else if (output_type == wall_times)
+ output_time = " wall time: " + wall.str() + ".";
+ else
+ output_time = ", CPU/wall time: " + cpu.str() + " / " + wall.str() + ".";
+
+ out_stream << section_names[section_to_exit] << output_time
+ << std::endl;
+ }
+
+ // delete the index from the list of
+ // active ones
+ std::vector<unsigned int>::iterator position_to_delete
+ = active_sections.begin() + active_index_to_delete;
+ active_sections.erase(position_to_delete);
+}
+
+
+
+void
+TimerOutput::print_summary ()
+{
+ // in case we want to write CPU times
+ if (output_type != wall_times)
+ {
+ double total_cpu_time;
+
+#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
+ {
+ double my_cpu_time = timer_all();
+ int mpiInitialized;
+ MPI_Initialized(&mpiInitialized);
+
+ if( mpiInitialized )
+ {
+ MPI_Allreduce (&my_cpu_time, &total_cpu_time, 1, MPI_DOUBLE, MPI_SUM,
+ MPI_COMM_WORLD);
+ }
+ else
+ total_cpu_time = my_cpu_time;
+ }
+#else
+ total_cpu_time = timer_all();
+#endif
+
+ // check that the sum of all times is
+ // less or equal than the total
+ // time. otherwise, we might have
+ // generated a lot of overhead in this
+ // function.
+ double check_time = 0.;
+ for (unsigned int i=0; i<section_names.size(); ++i)
+ {
+ check_time += section_total_cpu_times[i];
+ }
+ if (check_time > total_cpu_time)
+ {
+ total_cpu_time = check_time;
+ out_stream << std::endl << "Sum of counted times is larger than total time. "
+ << "Timer function may have introduced too much overhead." << std::endl;
+ }
+
+ // now generate a nice table
+ out_stream << "\n\n"
+ << "+---------------------------------------------+------------"
+ << "+------------+\n"
+ << "| Total CPU time elapsed since start |";
+ std::cout.width(10);
+ std::cout.precision(3);
+ out_stream << total_cpu_time << "s | |\n";
+ out_stream << "| | "
+ << "| |\n";
+ out_stream << "| Section | no. calls |";
+ std::cout.width(10);
+ std::cout.precision(3);
+ out_stream << " CPU time " << " | % of total |\n";
+ out_stream << "+---------------------------------------------+------------"
+ << "+------------+";
+ for (unsigned int i=0; i<section_names.size(); ++i)
+ {
+ std::string name_out = section_names[i];
+
+ // resize the array so that it is always
+ // of the same size
+ unsigned int pos_non_space = name_out.find_first_not_of (" ");
+ name_out.erase(0, pos_non_space);
+ name_out.resize (32, ' ');
+ out_stream << std::endl;
+ out_stream << "| " << name_out;
+ out_stream << "| ";
+ std::cout.width(9);
+ out_stream << section_n_calls[i] << " |";
+ std::cout.width(10);
+ std::cout.precision(3);
+ out_stream << section_total_cpu_times[i] << "s |";
+ std::cout.width(10);
+ std::cout.precision(2);
+ out_stream << section_total_cpu_times[i]/total_cpu_time * 100 << "% |";
+ }
+ out_stream << std::endl
+ << "+---------------------------------------------+"
+ << "------------+------------+\n"
+ << std::endl;
+ }
+
+ // in case we want to write out wallclock times
+ if (output_type != cpu_times)
+ {
+ double total_wall_time = timer_all.wall_time();
+
+ // check that the sum of all times is
+ // less or equal than the total
+ // time. otherwise, we might have
+ // generated a lot of overhead in this
+ // function.
+ double check_time = 0.;
+ for (unsigned int i=0; i<section_names.size(); ++i)
+ {
+ check_time += section_total_wall_times[i];
+ }
+ if (check_time > total_wall_time)
+ {
+ total_wall_time = check_time;
+ out_stream << std::endl
+ << "Sum of counted times is larger than total time. "
+ << "Timer function may have introduced too much overhead."
+ << std::endl;
+ }
+
+ // now generate a nice table
+ out_stream << "\n\n"
+ << "+---------------------------------------------+------------"
+ << "+------------+\n"
+ << "| Total wallclock time elapsed from start |";
+ std::cout.width(10);
+ std::cout.precision(3);
+ out_stream << total_wall_time << "s | |\n";
+ out_stream << "| | "
+ << "| |\n";
+ out_stream << "| Section | no. calls |";
+ std::cout.width(10);
+ std::cout.precision(3);
+ out_stream << " CPU time " << " | % of total |\n";
+ out_stream << "+---------------------------------------------+------------"
+ << "+------------+";
+ for (unsigned int i=0; i<section_names.size(); ++i)
+ {
+ std::string name_out = section_names[i];
+
+ // resize the array so that it is always
+ // of the same size
+ unsigned int pos_non_space = name_out.find_first_not_of (" ");
+ name_out.erase(0, pos_non_space);
+ name_out.resize (32, ' ');
+ out_stream << std::endl;
+ out_stream << "| " << name_out;
+ out_stream << "| ";
+ std::cout.width(9);
+ out_stream << section_n_calls[i] << " |";
+ std::cout.width(10);
+ std::cout.precision(3);
+ out_stream << section_total_wall_times[i] << "s |";
+ std::cout.width(10);
+ std::cout.precision(2);
+ out_stream << section_total_wall_times[i]/total_wall_time * 100 << "% |";
+ }
+ out_stream << std::endl
+ << "+---------------------------------------------+"
+ << "------------+------------+\n"
+ << std::endl;
+ }
+
+}
+
+
DEAL_II_NAMESPACE_CLOSE
#include <Epetra_Map.h>
+ // Time measurements.
+#include <base/timer.h>
+
#include <fstream>
#include <iostream>
#include <sstream>
// names into global namespace
using namespace dealii;
-
// @sect3{Equation data}
template <int dim>
double
- TemperatureInitialValues<dim>::value (const Point<dim> &p,
+ TemperatureInitialValues<dim>::value (const Point<dim> &,
const unsigned int) const
{
- return (p.norm() < 0.55+0.02*std::sin(p[0]*20) ? 1 : 0);
+ //return (p.norm() < 0.55+0.02*std::sin(p[0]*20) ? 1 : 0);
+ return 0.;
}
template <int dim>
double
- TemperatureRightHandSide<dim>::value (const Point<dim> &,
- const unsigned int /*component*/) const
+ TemperatureRightHandSide<dim>::value (const Point<dim> &p,
+ const unsigned int component) const
{
- return 0;
+ // return 0;
+ Assert (component == 0,
+ ExcMessage ("Invalid operation for a scalar function."));
+
+ Assert ((dim==2) || (dim==3), ExcNotImplemented());
+
+ static const Point<dim> source_centers[3]
+ = { (dim == 2 ? Point<dim>(.3,.1) : Point<dim>(.3,.5,.1)),
+ (dim == 2 ? Point<dim>(.45,.1) : Point<dim>(.45,.5,.1)),
+ (dim == 2 ? Point<dim>(.75,.1) : Point<dim>(.75,.5,.1)) };
+ static const double source_radius
+ = (dim == 2 ? 1./32 : 1./8);
+
+ return ((source_centers[0].distance (p) < source_radius)
+ ||
+ (source_centers[1].distance (p) < source_radius)
+ ||
+ (source_centers[2].distance (p) < source_radius)
+ ?
+ 1
+ :
+ 0);
+
}
void output_results () const;
void refine_mesh (const unsigned int max_grid_level);
- static
double
compute_viscosity(const std::vector<double> &old_temperature,
const std::vector<double> &old_old_temperature,
const std::vector<Tensor<1,dim> > &old_old_temperature_grads,
const std::vector<double> &old_temperature_laplacians,
const std::vector<double> &old_old_temperature_laplacians,
- const std::vector<Tensor<1,dim> > &present_velocity_values,
+ const std::vector<Tensor<1,dim> > &old_velocity_values,
+ const std::vector<Tensor<1,dim> > &old_old_velocity_values,
const std::vector<double> &gamma_values,
const double global_u_infty,
const double global_T_variation,
- const double global_Omega_diameter,
- const double cell_diameter,
- const double old_time_step);
+ const double cell_diameter);
const Epetra_Comm &trilinos_communicator;
ConditionalOStream pcout;
Triangulation<dim> triangulation;
+ double global_Omega_diameter;
const unsigned int stokes_degree;
FESystem<dim> stokes_fe;
TrilinosWrappers::BlockSparseMatrix stokes_preconditioner_matrix;
TrilinosWrappers::MPI::BlockVector stokes_solution;
+ TrilinosWrappers::BlockVector old_stokes_solution;
TrilinosWrappers::MPI::BlockVector stokes_rhs;
double old_time_step;
unsigned int timestep_number;
- boost::shared_ptr<TrilinosWrappers::PreconditionAMG> Amg_preconditioner;
- boost::shared_ptr<TrilinosWrappers::PreconditionIC> Mp_preconditioner;
+ boost::shared_ptr<TrilinosWrappers::PreconditionAMG> Amg_preconditioner;
+ boost::shared_ptr<TrilinosWrappers::PreconditionIC> Mp_preconditioner;
+ boost::shared_ptr<TrilinosWrappers::PreconditionIC> T_preconditioner;
bool rebuild_stokes_matrix;
- bool rebuild_temperature_matrices;
bool rebuild_stokes_preconditioner;
+ bool rebuild_temperature_matrices;
+ bool rebuild_temperature_preconditioner;
+
+ TimerOutput computing_timer;
};
old_time_step (0),
timestep_number (0),
rebuild_stokes_matrix (true),
+ rebuild_stokes_preconditioner (true),
rebuild_temperature_matrices (true),
- rebuild_stokes_preconditioner (true)
+ rebuild_temperature_preconditioner (true),
+ computing_timer (pcout, TimerOutput::summary,
+ TimerOutput::wall_times)
{}
template <int dim>
double
BoussinesqFlowProblem<dim>::
-compute_viscosity(const std::vector<double> &old_temperature,
- const std::vector<double> &old_old_temperature,
- const std::vector<Tensor<1,dim> > &old_temperature_grads,
- const std::vector<Tensor<1,dim> > &old_old_temperature_grads,
- const std::vector<double> &old_temperature_laplacians,
- const std::vector<double> &old_old_temperature_laplacians,
- const std::vector<Tensor<1,dim> > &present_velocity_values,
- const std::vector<double> &gamma_values,
- const double global_u_infty,
- const double global_T_variation,
- const double global_Omega_diameter,
- const double cell_diameter,
- const double old_time_step)
+compute_viscosity (const std::vector<double> &old_temperature,
+ const std::vector<double> &old_old_temperature,
+ const std::vector<Tensor<1,dim> > &old_temperature_grads,
+ const std::vector<Tensor<1,dim> > &old_old_temperature_grads,
+ const std::vector<double> &old_temperature_laplacians,
+ const std::vector<double> &old_old_temperature_laplacians,
+ const std::vector<Tensor<1,dim> > &old_velocity_values,
+ const std::vector<Tensor<1,dim> > &old_old_velocity_values,
+ const std::vector<double> &gamma_values,
+ const double global_u_infty,
+ const double global_T_variation,
+ const double cell_diameter)
{
- const double beta = 0.015 * dim;
- const double alpha = 1;
+ const double beta = 0.04 * dim;
+ const double alpha = 2;
if (global_u_infty == 0)
return 5e-3 * cell_diameter;
const unsigned int n_q_points = old_temperature.size();
- // Stage 1: calculate residual
double max_residual = 0;
double max_velocity = 0;
for (unsigned int q=0; q < n_q_points; ++q)
{
- const Tensor<1,dim> u = present_velocity_values[q];
+ const Tensor<1,dim> u = (old_velocity_values[q] +
+ old_old_velocity_values[q]) / 2;
const double dT_dt = (old_temperature[q] - old_old_temperature[q])
/ old_time_step;
return (beta *
max_velocity *
std::min (cell_diameter,
- std::pow(cell_diameter,alpha) * max_residual / global_scaling));
+ std::pow(cell_diameter,alpha) *
+ max_residual / global_scaling));
}
+
// @sect4{BoussinesqFlowProblem::setup_dofs}
template <int dim>
void BoussinesqFlowProblem<dim>::setup_dofs ()
{
+ computing_timer.enter_section("Setup dof systems");
std::vector<unsigned int> stokes_sub_blocks (dim+1,0);
stokes_sub_blocks[dim] = 1;
sp.compress();
stokes_matrix.reinit (sp);
+ /*std::cout << "Processor " << trilinos_communicator.MyPID()
+ << " stokes(0,0) rows: "
+ << stokes_matrix.block(0,0).matrix->NumMyRows()
+ << ", nnz: "
+ << stokes_matrix.block(0,0).matrix->NumMyNonzeros()
+ << std::endl;*/
+
}
{
0,
trilinos_communicator);
{
+ T_preconditioner.reset ();
temperature_mass_matrix.clear ();
temperature_stiffness_matrix.clear ();
temperature_matrix.clear ();
}
stokes_solution.reinit (stokes_partitioner);
+ old_stokes_solution.reinit (stokes_partitioner);
stokes_rhs.reinit (stokes_partitioner);
temperature_solution.reinit (temperature_partitioner);
old_temperature_solution.reinit (temperature_partitioner);
old_old_temperature_solution.reinit (temperature_partitioner);
temperature_rhs.reinit (temperature_partitioner);
+
+ computing_timer.exit_section();
}
if (rebuild_stokes_preconditioner == false)
return;
+ computing_timer.enter_section (" Build Stokes preconditioner");
pcout << " Rebuilding Stokes preconditioner..." << std::flush;
assemble_stokes_preconditioner ();
Amg_preconditioner->initialize(stokes_preconditioner_matrix.block(0,0),
TrilinosWrappers::PreconditionAMG::AdditionalData
- (true, true, 5e-2, null_space, 3, 0, false));
+ (true, true, 1e-2, null_space, 3, 0, false));
Mp_preconditioner = boost::shared_ptr<TrilinosWrappers::PreconditionIC>
(new TrilinosWrappers::PreconditionIC());
Mp_preconditioner->initialize (stokes_preconditioner_matrix.block(1,1));
- pcout << std::endl;
-
rebuild_stokes_preconditioner = false;
+
+ pcout << std::endl;
+ computing_timer.exit_section();
}
{
pcout << " Assembling..." << std::flush;
+ computing_timer.enter_section (" Assemble Stokes system");
+
if (rebuild_stokes_matrix == true)
stokes_matrix=0;
std::vector<SymmetricTensor<2,dim> > grads_phi_u (dofs_per_cell);
std::vector<double> div_phi_u (dofs_per_cell);
std::vector<double> phi_p (dofs_per_cell);
-
+
const FEValuesExtractors::Vector velocities (0);
const FEValuesExtractors::Scalar pressure (dim);
rebuild_stokes_matrix = false;
pcout << std::endl;
+ computing_timer.exit_section();
}
{
if (rebuild_temperature_matrices == false)
return;
-
+
+ computing_timer.enter_section (" Assemble temperature matrices");
temperature_mass_matrix = 0;
temperature_stiffness_matrix = 0;
temperature_stiffness_matrix.compress();
rebuild_temperature_matrices = false;
+ rebuild_temperature_preconditioner = true;
+
+ computing_timer.exit_section();
}
temperature_matrix.add (time_step, temperature_stiffness_matrix);
}
temperature_matrix.compress();
+
+ if (rebuild_temperature_preconditioner == true)
+ {
+ T_preconditioner = boost::shared_ptr<TrilinosWrappers::PreconditionIC>
+ (new TrilinosWrappers::PreconditionIC());
+ T_preconditioner->initialize (temperature_matrix);
+
+ rebuild_temperature_preconditioner = false;
+ }
temperature_rhs = 0;
std::vector<unsigned int> local_dof_indices (dofs_per_cell);
- std::vector<Tensor<1,dim> > present_velocity_values (n_q_points);
-
+ std::vector<Tensor<1,dim> > old_velocity_values (n_q_points);
+ std::vector<Tensor<1,dim> > old_old_velocity_values (n_q_points);
std::vector<double> old_temperature_values (n_q_points);
std::vector<double> old_old_temperature_values(n_q_points);
global_T_range = get_extrapolated_temperature_range();
const double global_Omega_diameter = GridTools::diameter (triangulation);
- const TrilinosWrappers::BlockVector localized_stokes_solution (stokes_solution);
+ const TrilinosWrappers::BlockVector
+ localized_stokes_solution (stokes_solution);
+ const TrilinosWrappers::BlockVector
+ localized_old_stokes_solution (old_stokes_solution);
const FEValuesExtractors::Vector velocities (0);
temperature_right_hand_side.value_list (temperature_fe_values.get_quadrature_points(),
gamma_values);
-
+
stokes_fe_values[velocities].get_function_values (localized_stokes_solution,
- present_velocity_values);
+ old_velocity_values);
+ stokes_fe_values[velocities].get_function_values (localized_old_stokes_solution,
+ old_old_velocity_values);
const double nu
= compute_viscosity (old_temperature_values,
old_old_temperature_grads,
old_temperature_laplacians,
old_old_temperature_laplacians,
- present_velocity_values,
+ old_velocity_values,
+ old_old_velocity_values,
gamma_values,
global_u_infty,
global_T_range.second - global_T_range.first,
- global_Omega_diameter, cell->diameter(),
- old_time_step);
+ cell->diameter());
for (unsigned int q=0; q<n_q_points; ++q)
{
phi_T[k] = temperature_fe_values.shape_value (k, q);
}
- const double old_T = old_temperature_values[q];
- const double old_old_T = old_old_temperature_values[q];
-
- const Tensor<1,dim> old_grad_T = old_temperature_grads[q];
- const Tensor<1,dim> old_old_grad_T = old_old_temperature_grads[q];
-
-
- const Tensor<1,dim> present_u = present_velocity_values[q];
-
- if (use_bdf2_scheme == true)
- {
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- local_rhs(i) += ((time_step + old_time_step) / old_time_step *
- old_T * phi_T[i]
- -
- (time_step * time_step) /
- (old_time_step * (time_step + old_time_step)) *
- old_old_T * phi_T[i]
- -
- time_step *
- present_u *
- ((1+time_step/old_time_step) * old_grad_T
- -
- time_step / old_time_step * old_old_grad_T) *
- phi_T[i]
- -
- time_step *
- nu *
- ((1+time_step/old_time_step) * old_grad_T
- -
- time_step / old_time_step * old_old_grad_T) *
- grad_phi_T[i]
- +
- time_step *
- gamma_values[q] * phi_T[i])
- *
- temperature_fe_values.JxW(q);
- }
- else
- {
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- local_rhs(i) += (old_T * phi_T[i]
- -
- time_step *
- present_u * old_grad_T * phi_T[i]
- -
- time_step *
- nu *
- old_grad_T * grad_phi_T[i]
- +
- time_step *
- gamma_values[q] * phi_T[i])
- *
- temperature_fe_values.JxW(q);
- }
+
+ const double old_Ts
+ = (use_bdf2_scheme ?
+ (old_temperature_values[q] *
+ (time_step + old_time_step) / old_time_step
+ -
+ old_old_temperature_values[q] *
+ (time_step * time_step) /
+ (old_time_step * (time_step + old_time_step)))
+ :
+ old_temperature_values[q]);
+
+ const Tensor<1,dim> ext_grad_T
+ = (use_bdf2_scheme ?
+ (old_temperature_grads[q] *
+ (1+time_step/old_time_step)
+ -
+ old_old_temperature_grads[q] *
+ time_step / old_time_step)
+ :
+ old_temperature_grads[q]);
+
+ const Tensor<1,dim> extrapolated_u
+ = (use_bdf2_scheme ?
+ (old_velocity_values[q] * (1+time_step/old_time_step) -
+ old_old_velocity_values[q] * time_step/old_time_step)
+ :
+ old_velocity_values[q]);
+
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ local_rhs(i) += (old_Ts * phi_T[i]
+ -
+ time_step *
+ extrapolated_u * ext_grad_T * phi_T[i]
+ -
+ time_step *
+ nu * ext_grad_T * grad_phi_T[i]
+ +
+ time_step *
+ gamma_values[q] * phi_T[i])
+ *
+ temperature_fe_values.JxW(q);
}
cell->get_dof_indices (local_dof_indices);
template <int dim>
void BoussinesqFlowProblem<dim>::solve ()
{
+ computing_timer.enter_section (" Solve Stokes system");
pcout << " Solving..." << std::endl;
{
stokes_constraints.distribute (localized_stokes_solution);
stokes_solution = localized_stokes_solution;
}
+ computing_timer.exit_section();
+
+
+ computing_timer.enter_section (" Assemble temperature rhs");
old_time_step = time_step;
- time_step = 1./(1.6*dim*std::sqrt(1.*dim)) /
+ time_step = 1./(1.9*dim*std::sqrt(1.*dim)) /
temperature_degree *
GridTools::minimal_cell_diameter(triangulation) /
- std::max (get_maximal_velocity(), 1.e-5);
+ std::max (get_maximal_velocity(), 0.01);
pcout << " " << "Time step: " << time_step
- << std::endl;
-
+ << std::endl;
+
temperature_solution = old_temperature_solution;
assemble_temperature_system ();
- {
+ computing_timer.exit_section ();
+
+ computing_timer.enter_section (" Solve temperature system");
+
+ {
SolverControl solver_control (temperature_matrix.m(),
1e-8*temperature_rhs.l2_norm());
SolverCG<TrilinosWrappers::MPI::Vector> cg (solver_control);
- TrilinosWrappers::PreconditionIC preconditioner;
- preconditioner.initialize (temperature_matrix);
-
cg.solve (temperature_matrix, temperature_solution,
- temperature_rhs, preconditioner);
+ temperature_rhs, *T_preconditioner);
TrilinosWrappers::Vector localized_temperature_solution (temperature_solution);
temperature_constraints.distribute (localized_temperature_solution);
pcout << " "
<< solver_control.last_step()
- << " CG iterations for temperature."
- << std::endl;
+ << " CG iterations for temperature" << std::endl;
+ computing_timer.exit_section();
double min_temperature = localized_temperature_solution(0),
max_temperature = localized_temperature_solution(0);
- for (unsigned int i=0; i<temperature_solution.size(); ++i)
+ for (unsigned int i=1; i<temperature_solution.size(); ++i)
{
min_temperature = std::min<double> (min_temperature,
localized_temperature_solution(i));
Assert (joint_dof_handler.n_dofs() ==
stokes_dof_handler.n_dofs() + temperature_dof_handler.n_dofs(),
ExcInternalError());
-
+
Vector<double> joint_solution (joint_dof_handler.n_dofs());
TrilinosWrappers::BlockVector localized_stokes_solution (stokes_solution);
TrilinosWrappers::Vector localized_temperature_solution (temperature_solution);
- {
- std::vector<unsigned int> local_joint_dof_indices (joint_fe.dofs_per_cell);
- std::vector<unsigned int> local_stokes_dof_indices (stokes_fe.dofs_per_cell);
- std::vector<unsigned int> local_temperature_dof_indices (temperature_fe.dofs_per_cell);
-
- typename DoFHandler<dim>::active_cell_iterator
- joint_cell = joint_dof_handler.begin_active(),
- joint_endc = joint_dof_handler.end(),
- stokes_cell = stokes_dof_handler.begin_active(),
- temperature_cell = temperature_dof_handler.begin_active();
- for (; joint_cell!=joint_endc; ++joint_cell, ++stokes_cell, ++temperature_cell)
+ if (Utilities::Trilinos::get_this_mpi_process(trilinos_communicator) == 0)
+ {
+
{
- joint_cell->get_dof_indices (local_joint_dof_indices);
- stokes_cell->get_dof_indices (local_stokes_dof_indices);
- temperature_cell->get_dof_indices (local_temperature_dof_indices);
-
- for (unsigned int i=0; i<joint_fe.dofs_per_cell; ++i)
- if (joint_fe.system_to_base_index(i).first.first == 0)
- {
- Assert (joint_fe.system_to_base_index(i).second
- <
- local_stokes_dof_indices.size(),
- ExcInternalError());
- joint_solution(local_joint_dof_indices[i])
- = localized_stokes_solution(local_stokes_dof_indices[joint_fe.system_to_base_index(i).second]);
- }
- else
- {
- Assert (joint_fe.system_to_base_index(i).first.first == 1,
- ExcInternalError());
- Assert (joint_fe.system_to_base_index(i).second
- <
- local_stokes_dof_indices.size(),
- ExcInternalError());
- joint_solution(local_joint_dof_indices[i])
- = localized_temperature_solution(local_temperature_dof_indices[joint_fe.system_to_base_index(i).second]);
- }
+ std::vector<unsigned int> local_joint_dof_indices (joint_fe.dofs_per_cell);
+ std::vector<unsigned int> local_stokes_dof_indices (stokes_fe.dofs_per_cell);
+ std::vector<unsigned int> local_temperature_dof_indices (temperature_fe.dofs_per_cell);
+
+ typename DoFHandler<dim>::active_cell_iterator
+ joint_cell = joint_dof_handler.begin_active(),
+ joint_endc = joint_dof_handler.end(),
+ stokes_cell = stokes_dof_handler.begin_active(),
+ temperature_cell = temperature_dof_handler.begin_active();
+ for (; joint_cell!=joint_endc; ++joint_cell, ++stokes_cell, ++temperature_cell)
+ {
+ joint_cell->get_dof_indices (local_joint_dof_indices);
+ stokes_cell->get_dof_indices (local_stokes_dof_indices);
+ temperature_cell->get_dof_indices (local_temperature_dof_indices);
+
+ for (unsigned int i=0; i<joint_fe.dofs_per_cell; ++i)
+ if (joint_fe.system_to_base_index(i).first.first == 0)
+ {
+ Assert (joint_fe.system_to_base_index(i).second
+ <
+ local_stokes_dof_indices.size(),
+ ExcInternalError());
+ joint_solution(local_joint_dof_indices[i])
+ = localized_stokes_solution(local_stokes_dof_indices[joint_fe.system_to_base_index(i).second]);
+ }
+ else
+ {
+ Assert (joint_fe.system_to_base_index(i).first.first == 1,
+ ExcInternalError());
+ Assert (joint_fe.system_to_base_index(i).second
+ <
+ local_stokes_dof_indices.size(),
+ ExcInternalError());
+ joint_solution(local_joint_dof_indices[i])
+ = localized_temperature_solution(local_temperature_dof_indices[joint_fe.system_to_base_index(i).second]);
+ }
+ }
}
- }
-
-
- std::vector<std::string> joint_solution_names (dim, "velocity");
- joint_solution_names.push_back ("p");
- joint_solution_names.push_back ("T");
- DataOut<dim> data_out;
+ std::vector<std::string> joint_solution_names (dim, "velocity");
+ joint_solution_names.push_back ("p");
+ joint_solution_names.push_back ("T");
+
+ DataOut<dim> data_out;
- if (Utilities::Trilinos::get_this_mpi_process(trilinos_communicator) == 0)
- {
data_out.attach_dof_handler (joint_dof_handler);
std::vector<DataComponentInterpretation::DataComponentInterpretation>
template <int dim>
void BoussinesqFlowProblem<dim>::refine_mesh (const unsigned int max_grid_level)
{
+ computing_timer.enter_section ("Refine mesh structure, part 1");
Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
TrilinosWrappers::Vector localized_temperature_solution (temperature_solution);
cell != triangulation.end(); ++cell)
cell->clear_refine_flag ();
- std::vector<TrilinosWrappers::Vector > x_solution (2);
- x_solution[0] = temperature_solution;
- x_solution[1] = old_temperature_solution;
+ std::vector<TrilinosWrappers::Vector> x_temperature (2);
+ x_temperature[0] = temperature_solution;
+ x_temperature[1] = old_temperature_solution;
+ TrilinosWrappers::BlockVector x_stokes = stokes_solution;
- SolutionTransfer<dim,TrilinosWrappers::Vector > soltrans(temperature_dof_handler);
+ SolutionTransfer<dim,TrilinosWrappers::Vector>
+ temperature_trans(temperature_dof_handler);
+ SolutionTransfer<dim,TrilinosWrappers::BlockVector>
+ stokes_trans(stokes_dof_handler);
triangulation.prepare_coarsening_and_refinement();
- soltrans.prepare_for_coarsening_and_refinement(x_solution);
+ temperature_trans.prepare_for_coarsening_and_refinement(x_temperature);
+ stokes_trans.prepare_for_coarsening_and_refinement(x_stokes);
triangulation.execute_coarsening_and_refinement ();
+ computing_timer.exit_section();
+
setup_dofs ();
- std::vector<TrilinosWrappers::Vector > tmp (2);
- tmp[0] = temperature_solution;
- tmp[1] = temperature_solution;
- soltrans.interpolate(x_solution, tmp);
+ computing_timer.enter_section ("Refine mesh structure, part 2");
+
+ std::vector<TrilinosWrappers::Vector> tmp (2);
+ tmp[0].reinit (temperature_solution);
+ tmp[1].reinit (temperature_solution);
+ temperature_trans.interpolate(x_temperature, tmp);
temperature_solution = tmp[0];
old_temperature_solution = tmp[1];
- rebuild_stokes_matrix = true;
- rebuild_temperature_matrices = true;
- rebuild_stokes_preconditioner = true;
+ TrilinosWrappers::BlockVector x_stokes_new = stokes_solution;
+ stokes_trans.interpolate (x_stokes, x_stokes_new);
+ stokes_solution = x_stokes_new;
+
+ rebuild_stokes_matrix = true;
+ rebuild_stokes_preconditioner = true;
+ rebuild_temperature_matrices = true;
+ rebuild_temperature_preconditioner = true;
+
+ computing_timer.exit_section();
}
const unsigned int initial_refinement = (dim == 2 ? 3 : 2);
const unsigned int n_pre_refinement_steps = (dim == 2 ? 4 : 3);
- GridGenerator::half_hyper_shell (triangulation,
- Point<dim>(), 0.5, 1.0);
+ //GridGenerator::half_hyper_shell (triangulation,
+ // Point<dim>(), 0.5, 1.0);
- static HalfHyperShellBoundary<dim> boundary;
- triangulation.set_boundary (0, boundary);
+ //static HyperShellBoundary<dim> boundary;
+ //triangulation.set_boundary (0, boundary);
+ GridGenerator::hyper_cube (triangulation);
+ global_Omega_diameter = GridTools::diameter (triangulation);
triangulation.refine_global (initial_refinement);
start_time_iteration:
- VectorTools::project (temperature_dof_handler,
- temperature_constraints,
- QGauss<dim>(temperature_degree+2),
- EquationData::TemperatureInitialValues<dim>(),
- old_temperature_solution);
+ //VectorTools::project (temperature_dof_handler,
+ // temperature_constraints,
+ // QGauss<dim>(temperature_degree+2),
+ // EquationData::TemperatureInitialValues<dim>(),
+ // old_temperature_solution);
+ assemble_temperature_matrix ();
+ // Create right hand side for projection
+ {
+ QGauss<dim> quadrature(temperature_degree+2);
+ UpdateFlags update_flags = UpdateFlags(update_values |
+ update_quadrature_points |
+ update_JxW_values);
+ FEValues<dim> fe_values (temperature_fe, quadrature, update_flags);
+
+ const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
+ n_q_points = fe_values.n_quadrature_points;
+
+ std::vector<unsigned int> dofs (dofs_per_cell);
+ Vector<double> cell_vector (dofs_per_cell);
+
+ typename DoFHandler<dim>::active_cell_iterator
+ cell = temperature_dof_handler.begin_active(),
+ endc = temperature_dof_handler.end();
+
+ std::vector<double> rhs_values(n_q_points);
+
+ TrilinosWrappers::MPI::Vector rhs (temperature_mass_matrix.matrix->RowMap()),
+ sol (temperature_mass_matrix.matrix->RowMap());
+
+ for (; cell!=endc; ++cell)
+ if (cell->subdomain_id() ==
+ Utilities::Trilinos::get_this_mpi_process(trilinos_communicator))
+ {
+ fe_values.reinit(cell);
+
+ const std::vector<double> &weights = fe_values.get_JxW_values ();
+ EquationData::TemperatureInitialValues<dim>().value_list
+ (fe_values.get_quadrature_points(), rhs_values);
+
+ cell_vector = 0;
+ for (unsigned int point=0; point<n_q_points; ++point)
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ cell_vector(i) += rhs_values[point] *
+ fe_values.shape_value(i,point) *
+ weights[point];
+
+ cell->get_dof_indices (dofs);
+
+ temperature_constraints.distribute_local_to_global (cell_vector,
+ dofs,
+ rhs);
+ }
+
+ ReductionControl control(5*rhs.size(), 0., 1e-12, false, false);
+ GrowingVectorMemory<TrilinosWrappers::MPI::Vector> memory;
+ SolverCG<TrilinosWrappers::MPI::Vector> cg(control,memory);
+
+ TrilinosWrappers::PreconditionIC prec;
+ prec.initialize(temperature_mass_matrix);
+ // solve
+ cg.solve (temperature_mass_matrix, sol, rhs, prec);
+
+ // distribute solution
+ old_temperature_solution = sol;
+ temperature_constraints.distribute (old_temperature_solution);
+ }
timestep_number = 0;
time_step = old_time_step = 0;
time += time_step;
++timestep_number;
+ old_stokes_solution = stokes_solution;
old_old_temperature_solution = old_temperature_solution;
old_temperature_solution = temperature_solution;
}
- while (time <= 100);
+ while (time <= 0);
+
}
Utilities::System::MPI_InitFinalize mpi_initialization(argc, argv);
- BoussinesqFlowProblem<2> flow_problem;
+ BoussinesqFlowProblem<3> flow_problem;
flow_problem.run ();
}
catch (std::exception &exc)
projects / dealii-svn.git / commitdiff