[cig-commits] commit 1863 by buerg to /var/svn/dealii/aspect
dealii.demon at gmail.com
dealii.demon at gmail.com
Tue Aug 27 08:34:28 PDT 2013
Revision 1863
Assemble composition systems at once.
U trunk/aspire/include/aspect/simulator.h
U trunk/aspire/source/simulator/assembly.cc
U trunk/aspire/source/simulator/core.cc
http://www.dealii.org/websvn/revision.php?repname=Aspect+Repository&path=%2F&rev=1863&peg=1863
Diff:
Modified: trunk/aspire/include/aspect/simulator.h
===================================================================
--- trunk/aspire/include/aspect/simulator.h 2013-08-27 15:24:02 UTC (rev 1862)
+++ trunk/aspire/include/aspect/simulator.h 2013-08-27 15:33:15 UTC (rev 1863)
@@ -423,7 +423,7 @@
* This function is implemented in
* <code>source/simulator/assembly.cc</code>.
*/
- void assemble_composition_system (const unsigned int c);
+ void assemble_composition_systems ();
/**
* Solve the temperature linear system. Return the initial nonlinear residual, i.e.,
@@ -620,15 +620,14 @@
* <code>source/simulator/assembly.cc</code>.
*/
void
- local_assemble_composition_system (const unsigned int c,
- const std::pair<double,double> global_field_range,
- const double global_max_velocity,
- const double global_entropy_variation,
- const SynchronousIterators<std_cxx1x::tuple<FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>,
- FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>,
- FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> > > &cell,
- internal::Assembly::Scratch::CompositionSystem<dim> &scratch,
- internal::Assembly::CopyData::CompositionSystem<dim> &data);
+ local_assemble_composition_systems (const std::vector<std::pair<double,double> > global_field_ranges,
+ const double global_max_velocity,
+ const std::vector<double> global_entropy_variations,
+ const SynchronousIterators<std_cxx1x::tuple<FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>,
+ FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>,
+ FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> > > &cell,
+ internal::Assembly::Scratch::CompositionSystem<dim> &scratch,
+ internal::Assembly::CopyData::CompositionSystem<dim> &data);
/**
@@ -650,8 +649,7 @@
* <code>source/simulator/assembly.cc</code>.
*/
void
- copy_local_to_global_composition_system (const unsigned int c,
- const internal::Assembly::CopyData::CompositionSystem<dim> &data);
+ copy_local_to_global_composition_systems (const internal::Assembly::CopyData::CompositionSystem<dim> &data);
/**
* @}
@@ -911,8 +909,6 @@
std::vector<ConstraintMatrix> current_constraints_composition;
double pressure_adjustment;
- double pressure_scaling;
- double pressure_scaling_old;
/**
* @}
Modified: trunk/aspire/source/simulator/assembly.cc
===================================================================
--- trunk/aspire/source/simulator/assembly.cc 2013-08-27 15:24:02 UTC (rev 1862)
+++ trunk/aspire/source/simulator/assembly.cc 2013-08-27 15:33:15 UTC (rev 1863)
@@ -24,6 +24,7 @@
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/work_stream.h>
+#include <deal.II/lac/block_vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/grid/tria_iterator.h>
@@ -193,17 +194,17 @@
const unsigned int n_compositional_fields);
CompositionSystem (const CompositionSystem &data);
- FEValues<dim> fe_values_velocity;
- FEValues<dim> fe_values_temperature;
- FEValues<dim> fe_values_composition;
+ FEValues<dim> fe_values_velocity;
+ FEValues<dim> fe_values_temperature;
+ FEValues<dim> fe_values_composition;
- std::vector<double> phi_field;
- std::vector<Tensor<1,dim> > grad_phi_field;
+ std::vector<double> phi_field;
+ std::vector<Tensor<1,dim> > grad_phi_field;
- std::vector<double> old_composition_values;
- std::vector<double> old_old_composition_values;
+ std::vector<std::vector<double> > old_composition_values;
+ std::vector<double> old_old_composition_values;
- std::vector<Tensor<1,dim> > current_velocity_values;
+ std::vector<Tensor<1,dim> > current_velocity_values;
std::vector<std::vector<double> > current_composition_values;
typename MaterialModel::Interface<dim>::MaterialModelInputs material_model_inputs;
@@ -272,7 +273,7 @@
phi_field (fe_composition.dofs_per_cell),
grad_phi_field (fe_composition.dofs_per_cell),
- old_composition_values(quadrature.size ()),
+ old_composition_values (n_compositional_fields, std::vector<double> (quadrature.size ())),
old_old_composition_values(quadrature.size ()),
current_velocity_values(quadrature.size()),
current_composition_values(n_compositional_fields,
@@ -334,7 +335,7 @@
phi_field (scratch.phi_field),
grad_phi_field (scratch.grad_phi_field),
- old_composition_values(scratch.old_composition_values),
+ old_composition_values (scratch.old_composition_values),
old_old_composition_values(scratch.old_old_composition_values),
current_velocity_values(scratch.current_velocity_values),
current_composition_values(scratch.current_composition_values),
@@ -439,11 +440,11 @@
template <int dim>
struct CompositionSystem
{
- CompositionSystem (const FiniteElement<dim> &finite_element);
+ CompositionSystem (const FiniteElement<dim> &finite_element, const unsigned int n_compositional_fields);
CompositionSystem (const CompositionSystem &data);
- FullMatrix<double> local_matrix;
- Vector<double> local_rhs;
+ std::vector<FullMatrix<double> > local_matrices;
+ BlockVector<double> local_rhs;
std::vector<types::global_dof_index> local_dof_indices;
};
@@ -463,11 +464,11 @@
template <int dim>
CompositionSystem<dim>::
- CompositionSystem (const FiniteElement<dim> &finite_element)
+ CompositionSystem (const FiniteElement<dim> &finite_element, const unsigned int n_compositional_fields)
:
- local_matrix (finite_element.dofs_per_cell,
- finite_element.dofs_per_cell),
- local_rhs (finite_element.dofs_per_cell),
+ local_matrices (n_compositional_fields, FullMatrix<double> (finite_element.dofs_per_cell,
+ finite_element.dofs_per_cell)),
+ local_rhs (n_compositional_fields, finite_element.dofs_per_cell),
local_dof_indices (finite_element.dofs_per_cell)
{}
@@ -487,7 +488,7 @@
CompositionSystem<dim>::
CompositionSystem (const CompositionSystem &data)
:
- local_matrix (data.local_matrix),
+ local_matrices (data.local_matrices),
local_rhs (data.local_rhs),
local_dof_indices (data.local_dof_indices)
{}
@@ -818,8 +819,7 @@
// for (unsigned int n = 0; n < parameters.n_compositional_fields; ++n)
// right_hand_side += (scratch.velocity_values[q] * scratch.composition_gradients[n][q]) / scratch.composition_values[n][q];
- right_hand_side *= pressure_scaling * scratch.fe_values_velocity.JxW (q)
- * (axisymmetric_terms ? scratch.fe_values_velocity.quadrature_point (q) (0) : 1);
+ right_hand_side *= scratch.fe_values_velocity.JxW (q) * (axisymmetric_terms ? scratch.fe_values_velocity.quadrature_point (q) (0) : 1);
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
@@ -830,10 +830,9 @@
? -2.0/3.0 * eta + parameters.stabilization_graddiv
: parameters.stabilization_graddiv)
* (scratch.div_phi_u[i] * scratch.div_phi_u[j])
- + pressure_scaling *
- (pressure_scaling * scratch.phi_p[i] * scratch.phi_p[j]
- - scratch.div_phi_u[i] * scratch.phi_p[j]
- - scratch.phi_p[i] * scratch.div_phi_u[j]))
+ + scratch.phi_p[i] * scratch.phi_p[j]
+ - scratch.div_phi_u[i] * scratch.phi_p[j]
+ - scratch.phi_p[i] * scratch.div_phi_u[j])
*
(axisymmetric_terms ?
scratch.fe_values_velocity.quadrature_point (q) (0) : 1)
@@ -1053,25 +1052,29 @@
template <int dim>
void Simulator<dim>::
- local_assemble_composition_system (const unsigned int c,
- const std::pair<double,double> global_field_range,
- const double global_max_velocity,
- const double global_entropy_variation,
- const SynchronousIterators<std_cxx1x::tuple<FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>,
- FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>,
- FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> > > &cell,
- internal::Assembly::Scratch::CompositionSystem<dim> &scratch,
- internal::Assembly::CopyData::CompositionSystem<dim> &data)
+ local_assemble_composition_systems (const std::vector<std::pair<double, double> > global_field_ranges,
+ const double global_max_velocity,
+ const std::vector<double> global_entropy_variations,
+ const SynchronousIterators<std_cxx1x::tuple<FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>,
+ FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>,
+ FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> > > &cell,
+ internal::Assembly::Scratch::CompositionSystem<dim> &scratch,
+ internal::Assembly::CopyData::CompositionSystem<dim> &data)
{
const unsigned int dofs_per_cell = scratch.fe_values_composition.get_fe ().dofs_per_cell;
const unsigned int n_q_points = scratch.fe_values_composition.n_quadrature_points;
scratch.fe_values_composition.reinit (std_cxx1x::get<2> (cell.iterators));
std_cxx1x::get<2> (cell.iterators)->get_dof_indices (data.local_dof_indices);
- data.local_matrix = 0;
+
+ for (unsigned int c = 0; c < parameters.n_compositional_fields; ++c)
+ {
+ data.local_matrices[c] = 0;
+ scratch.fe_values_composition.get_function_values (old_solution_composition.block (c),
+ scratch.old_composition_values[c]);
+ }
+
data.local_rhs = 0;
- scratch.fe_values_composition.get_function_values (old_solution_composition.block (c),
- scratch.old_composition_values);
scratch.fe_values_velocity.reinit (std_cxx1x::get<0> (cell.iterators));
scratch.fe_values_velocity[introspection.extractors.velocities].get_function_values (current_linearization_point_velocity,
scratch.current_velocity_values);
@@ -1094,29 +1097,33 @@
scratch.grad_phi_field[i] = scratch.fe_values_composition.shape_grad (i, q);
scratch.phi_field[i] = scratch.fe_values_composition.shape_value (i, q);
}
-
- const double old_value = scratch.old_composition_values[q];
+
const double rho = scratch.material_model_outputs.densities[q];
- const double diffusivity = rho * scratch.material_model_outputs.diffusivities[q];
- const double source_term = scratch.material_model_outputs.compositional_sources[q][c];
const Tensor<1, dim> current_u = scratch.current_velocity_values[q];
+
+ for (unsigned int c = 0; c < parameters.n_compositional_fields; ++c)
+ {
+ const double diffusivity = rho * scratch.material_model_outputs.diffusivities[q];
+ const double old_value = scratch.old_composition_values[c][q];
+ const double source_term = scratch.material_model_outputs.compositional_sources[q][c];
- for (unsigned int i = 0; i < dofs_per_cell; ++i)
- {
- for (unsigned int j = 0; j < dofs_per_cell; ++j)
- {
- data.local_matrix (i, j)
- += (time_step * (diffusivity * (scratch.grad_phi_field[i] * scratch.grad_phi_field[j])
- + rho * scratch.phi_field[i] * (current_u * scratch.grad_phi_field[j]))
- + scratch.phi_field[i] * scratch.phi_field[j])
- * scratch.fe_values_composition.JxW (q)
- * (axisymmetric_terms ?
- scratch.fe_values_composition.quadrature_point (q) (0) : 1.0);
- }
+ for (unsigned int i = 0; i < dofs_per_cell; ++i)
+ {
+ for (unsigned int j = 0; j < dofs_per_cell; ++j)
+ {
+ data.local_matrices[c] (i, j)
+ += (time_step * (diffusivity * (scratch.grad_phi_field[i] * scratch.grad_phi_field[j])
+ + rho * scratch.phi_field[i] * (current_u * scratch.grad_phi_field[j]))
+ + scratch.phi_field[i] * scratch.phi_field[j])
+ * scratch.fe_values_composition.JxW (q)
+ * (axisymmetric_terms ?
+ scratch.fe_values_composition.quadrature_point (q) (0) : 1.0);
+ }
- data.local_rhs (i) += scratch.phi_field[i] * (time_step * source_term + old_value) * scratch.fe_values_composition.JxW (q)
- * (axisymmetric_terms ? scratch.fe_values_composition.quadrature_point (q) (0) : 1.0);
- }
+ data.local_rhs.block (c) (i) += scratch.phi_field[i] * (time_step * source_term + old_value) * scratch.fe_values_composition.JxW (q)
+ * (axisymmetric_terms ? scratch.fe_values_composition.quadrature_point (q) (0) : 1.0);
+ }
+ }
}
}
@@ -1137,14 +1144,14 @@
template <int dim>
void
Simulator<dim>::
- copy_local_to_global_composition_system (const unsigned int c,
- const internal::Assembly::CopyData::CompositionSystem<dim> &data)
+ copy_local_to_global_composition_systems (const internal::Assembly::CopyData::CompositionSystem<dim> &data)
{
- current_constraints_composition[c].distribute_local_to_global (data.local_matrix,
- data.local_rhs,
- data.local_dof_indices,
- system_matrix_composition.block (c, c),
- system_rhs_composition.block (c));
+ for (unsigned int c = 0; c < parameters.n_compositional_fields; ++c)
+ current_constraints_composition[c].distribute_local_to_global (data.local_matrices[c],
+ data.local_rhs.block (c),
+ data.local_dof_indices,
+ system_matrix_composition.block (c, c),
+ system_rhs_composition.block (c));
}
@@ -1218,13 +1225,20 @@
template <int dim>
- void Simulator<dim>::assemble_composition_system (const unsigned int c)
+ void Simulator<dim>::assemble_composition_systems ()
{
computing_timer.enter_section (" Assemble composition system");
- system_matrix_composition.block (c, c) = 0.0;
- system_rhs_composition.block (c) = 0.0;
- const std::pair<double, double> global_field_range = get_extrapolated_composition_range (c);
+ std::vector<double> global_entropy_variations (parameters.n_compositional_fields);
+ std::vector<std::pair<double, double> > global_field_ranges (parameters.n_compositional_fields);
+
+ for (unsigned int c = 0; c < parameters.n_compositional_fields; ++c)
+ {
+ global_field_ranges[c] = get_extrapolated_composition_range (c);
+ global_entropy_variations[c] = get_composition_variation (0.5 * (global_field_ranges[c].first + global_field_ranges[c].second), c);
+ system_matrix_composition.block (c, c) = 0.0;
+ system_rhs_composition.block (c) = 0.0;
+ }
FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
cell_filter_velocity_begin (IteratorFilters::LocallyOwnedCell (), dof_handler_velocity.begin_active ());
@@ -1255,34 +1269,33 @@
run (filtered_iterators_begin,
filtered_iterators_end,
std_cxx1x::bind (&Simulator<dim>::
- local_assemble_composition_system,
+ local_assemble_composition_systems,
this,
- c,
- global_field_range,
+ global_field_ranges,
get_maximal_velocity (old_solution_velocity),
// use the mid-value of the advected field instead of the
// integral mean. results are not very
// sensitive to this and this is far simpler
- get_composition_variation (0.5 * (global_field_range.first +
- global_field_range.second),
- c),
+ global_entropy_variations,
std_cxx1x::_1,
std_cxx1x::_2,
std_cxx1x::_3),
std_cxx1x::bind (&Simulator<dim>::
- copy_local_to_global_composition_system,
+ copy_local_to_global_composition_systems,
this,
- c,
std_cxx1x::_1),
internal::Assembly::Scratch::
CompositionSystem<dim> (fe_velocity, fe_temperature, fe_composition, mapping,
QGauss<dim>(parameters.composition_degree+2),
parameters.n_compositional_fields),
internal::Assembly::CopyData::
- CompositionSystem<dim> (fe_composition));
+ CompositionSystem<dim> (fe_composition, parameters.n_compositional_fields));
- system_matrix_composition.block (c, c).compress (VectorOperation::add);
- system_rhs_composition.block (c).compress (VectorOperation::add);
+ for (unsigned int c = 0; c < parameters.n_compositional_fields; ++c)
+ {
+ system_matrix_composition.block (c, c).compress (VectorOperation::add);
+ system_rhs_composition.block (c).compress (VectorOperation::add);
+ }
computing_timer.exit_section();
}
@@ -1313,20 +1326,18 @@
FilteredIterator<DoFHandler<dim>::active_cell_iterator> > > &cell, \
internal::Assembly::Scratch::TemperatureSystem<dim> &scratch, \
internal::Assembly::CopyData::TemperatureSystem<dim> &data); \
- template void Simulator<dim>::local_assemble_composition_system (const unsigned int c, \
- const std::pair<double,double> global_field_range, \
- const double global_max_velocity, \
- const double global_entropy_variation, \
- const SynchronousIterators<std_cxx1x::tuple<FilteredIterator<DoFHandler<dim>::active_cell_iterator>, \
- FilteredIterator<DoFHandler<dim>::active_cell_iterator>, \
- FilteredIterator<DoFHandler<dim>::active_cell_iterator> > > &cell, \
- internal::Assembly::Scratch::CompositionSystem<dim> &scratch, \
- internal::Assembly::CopyData::CompositionSystem<dim> &data); \
+ template void Simulator<dim>::local_assemble_composition_systems (const std::vector<std::pair<double, double> > global_field_ranges, \
+ const double global_max_velocity, \
+ const std::vector<double> global_entropy_variations, \
+ const SynchronousIterators<std_cxx1x::tuple<FilteredIterator<DoFHandler<dim>::active_cell_iterator>, \
+ FilteredIterator<DoFHandler<dim>::active_cell_iterator>, \
+ FilteredIterator<DoFHandler<dim>::active_cell_iterator> > > &cell, \
+ internal::Assembly::Scratch::CompositionSystem<dim> &scratch, \
+ internal::Assembly::CopyData::CompositionSystem<dim> &data); \
template void Simulator<dim>::copy_local_to_global_temperature_system (const internal::Assembly::CopyData::TemperatureSystem<dim> &data); \
- template void Simulator<dim>::copy_local_to_global_composition_system (const unsigned int c, \
- const internal::Assembly::CopyData::CompositionSystem<dim> &data); \
+ template void Simulator<dim>::copy_local_to_global_composition_systems (const internal::Assembly::CopyData::CompositionSystem<dim> &data); \
template void Simulator<dim>::assemble_temperature_system (); \
- template void Simulator<dim>::assemble_composition_system (const unsigned int c); \
+ template void Simulator<dim>::assemble_composition_systems (); \
Modified: trunk/aspire/source/simulator/core.cc
===================================================================
--- trunk/aspire/source/simulator/core.cc 2013-08-27 15:24:02 UTC (rev 1862)
+++ trunk/aspire/source/simulator/core.cc 2013-08-27 15:33:15 UTC (rev 1863)
@@ -241,8 +241,6 @@
*geometry_model);
compositional_boundary_conditions[p->first].reset (bv);
}
-
- pressure_scaling = material_model->reference_viscosity() / geometry_model->length_scale();
// make sure that we don't have to fill every column of the statistics
// object in each time step.
@@ -1041,10 +1039,10 @@
build_temperature_preconditioner (T_preconditioner);
solve_temperature ();
current_linearization_point_temperature = solution_temperature;
+ assemble_composition_systems ();
for (unsigned int c = 0; c < parameters.n_compositional_fields; ++c)
{
- assemble_composition_system (c);
build_composition_preconditioner (c, C_preconditioner);
solve_composition (c);
}
More information about the CIG-COMMITS
mailing list