[cig-commits] commit 1948 by bangerth to /var/svn/dealii/aspect
dealii.demon at gmail.com
dealii.demon at gmail.com
Tue Oct 8 15:38:18 PDT 2013
Revision 1948
Only build that submatrix in advection assembly that corresponds to the actual advection variables we care about. Compress a couple of arrays in the scratch data structure that then no longer need to be so large (currently 81+8+27=116 elements, after this patch only 27 elements).
U trunk/aspect/source/simulator/assembly.cc
http://www.dealii.org/websvn/revision.php?repname=Aspect+Repository&path=%2F&rev=1948&peg=1948
Diff:
Modified: trunk/aspect/source/simulator/assembly.cc
===================================================================
--- trunk/aspect/source/simulator/assembly.cc 2013-10-08 21:49:10 UTC (rev 1947)
+++ trunk/aspect/source/simulator/assembly.cc 2013-10-08 22:37:29 UTC (rev 1948)
@@ -181,6 +181,7 @@
struct AdvectionSystem
{
AdvectionSystem (const FiniteElement<dim> &finite_element,
+ const FiniteElement<dim> &advection_element,
const Mapping<dim> &mapping,
const Quadrature<dim> &quadrature,
const unsigned int n_compositional_fields);
@@ -188,6 +189,16 @@
FEValues<dim> finite_element_values;
+ std::vector<types::global_dof_index> local_dof_indices;
+
+ /**
+ * Variables describing the values and gradients of the
+ * shape functions at the quadrature points, as they are
+ * used in the advection assembly function. note that the sizes
+ * of these arrays are equal to the number of shape functions
+ * corresponding to the advected field (and not of the overall
+ * field!), and that they are also correspondingly indexed.
+ */
std::vector<double> phi_field;
std::vector<Tensor<1,dim> > grad_phi_field;
@@ -231,6 +242,7 @@
template <int dim>
AdvectionSystem<dim>::
AdvectionSystem (const FiniteElement<dim> &finite_element,
+ const FiniteElement<dim> &advection_element,
const Mapping<dim> &mapping,
const Quadrature<dim> &quadrature,
const unsigned int n_compositional_fields)
@@ -243,8 +255,10 @@
update_quadrature_points |
update_JxW_values),
- phi_field (finite_element.dofs_per_cell),
- grad_phi_field (finite_element.dofs_per_cell),
+ local_dof_indices (finite_element.dofs_per_cell),
+
+ phi_field (advection_element.dofs_per_cell),
+ grad_phi_field (advection_element.dofs_per_cell),
old_velocity_values (quadrature.size()),
old_old_velocity_values (quadrature.size()),
old_pressure (quadrature.size()),
@@ -284,6 +298,8 @@
scratch.finite_element_values.get_quadrature(),
scratch.finite_element_values.get_update_flags()),
+ local_dof_indices (scratch.finite_element_values.get_fe().dofs_per_cell),
+
phi_field (scratch.phi_field),
grad_phi_field (scratch.grad_phi_field),
old_velocity_values (scratch.old_velocity_values),
@@ -393,11 +409,31 @@
template <int dim>
struct AdvectionSystem
{
+ /**
+ * Constructor.
+ * @param finite_element The element that describes the field for which we
+ * are trying to assemble a linear system. <b>Not</b> the global finite
+ * element.
+ */
AdvectionSystem (const FiniteElement<dim> &finite_element);
AdvectionSystem (const AdvectionSystem &data);
+ /**
+ * Local contributions to the global matrix and right hand side
+ * that correspond only to the variables listed in local_dof_indices
+ */
FullMatrix<double> local_matrix;
Vector<double> local_rhs;
+
+ /**
+ * Indices of those degrees of freedom that actually correspond
+ * to the temperature or compositional field. since this structure
+ * is used to represent just contributions to the advection
+ * systems, there will be no contributions to other parts of the
+ * system and consequently, we do not need to list here indices
+ * that correspond to velocity or pressure degrees (or, in fact
+ * any other variable outside the block we are currently considering)
+ */
std::vector<types::global_dof_index> local_dof_indices;
};
@@ -569,7 +605,11 @@
const double field = ((*scratch.old_field_values)[q] + (*scratch.old_old_field_values)[q]) / 2;
const double gamma
- = compute_heating_term(scratch,scratch.explicit_material_model_inputs,scratch.explicit_material_model_outputs,temperature_or_composition,q);
+ = compute_heating_term(scratch,
+ scratch.explicit_material_model_inputs,
+ scratch.explicit_material_model_outputs,
+ temperature_or_composition,
+ q);
double residual
= std::abs(density * c_P * (dField_dt + u_grad_field) - k_Delta_field - gamma);
@@ -1149,12 +1189,19 @@
const typename DoFHandler<dim>::active_cell_iterator &cell,
internal::Assembly::Scratch::AdvectionSystem<dim> &scratch,
internal::Assembly::CopyData::AdvectionSystem<dim> &data)
-{
+ {
const bool use_bdf2_scheme = (timestep_number > 1);
- const unsigned int dofs_per_cell = scratch.finite_element_values.get_fe().dofs_per_cell;
const unsigned int n_q_points = scratch.finite_element_values.n_quadrature_points;
+ // also have the number of dofs that correspond just to the element for
+ // the system we are currently trying to assemble
+ const unsigned int advection_dofs_per_cell = data.local_dof_indices.size();
+
+ Assert (advection_dofs_per_cell < scratch.finite_element_values.get_fe().dofs_per_cell, ExcInternalError());
+ Assert (scratch.grad_phi_field.size() == advection_dofs_per_cell, ExcInternalError());
+ Assert (scratch.phi_field.size() == advection_dofs_per_cell, ExcInternalError());
+
const unsigned int solution_component
= (temperature_or_composition.is_temperature()
?
@@ -1172,8 +1219,13 @@
);
scratch.finite_element_values.reinit (cell);
- cell->get_dof_indices (data.local_dof_indices);
+ // get all dof indices on the current cell, then extract those
+ // that correspond to the solution_field we are interested in
+ cell->get_dof_indices (scratch.local_dof_indices);
+ for (unsigned int i=0; i<advection_dofs_per_cell; ++i)
+ data.local_dof_indices[i] = scratch.local_dof_indices[scratch.finite_element_values.get_fe().component_to_system_index(solution_component, i)];
+
data.local_matrix = 0;
data.local_rhs = 0;
@@ -1265,14 +1317,14 @@
for (unsigned int q=0; q<n_q_points; ++q)
{
- for (unsigned int k=0; k<dofs_per_cell; ++k)
- // We only need to look up values of shape functions if they
- // belong to 'our' component. They are zero otherwise anyway.
- // Note that we later only look at the values that we do set here.
- if (cell->get_fe().system_to_component_index(k).first == solution_component)
+ // precompute the values of shape functions and their gradients.
+ // We only need to look up values of shape functions if they
+ // belong to 'our' component. They are zero otherwise anyway.
+ // Note that we later only look at the values that we do set here.
+ for (unsigned int k=0; k<advection_dofs_per_cell; ++k)
{
- scratch.grad_phi_field[k] = scratch.finite_element_values[solution_field].gradient (k,q);
- scratch.phi_field[k] = scratch.finite_element_values[solution_field].value (k, q);
+ scratch.grad_phi_field[k] = scratch.finite_element_values[solution_field].gradient (scratch.finite_element_values.get_fe().component_to_system_index(solution_component, k),q);
+ scratch.phi_field[k] = scratch.finite_element_values[solution_field].value (scratch.finite_element_values.get_fe().component_to_system_index(solution_component, k), q);
}
const double density_c_P =
@@ -1312,18 +1364,19 @@
const double factor = (use_bdf2_scheme)? ((2*time_step + old_time_step) /
(time_step + old_time_step)) : 1.0;
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- if (cell->get_fe().system_to_component_index(i).first == solution_component)
+ // do the actual assembly. note that we only need to loop over the advection
+ // shape functions because these are the only contributions we compute here
+ for (unsigned int i=0; i<advection_dofs_per_cell; ++i)
{
- data.local_rhs(i) += (field_term_for_rhs * scratch.phi_field[i]
+ data.local_rhs(i)
+ += (field_term_for_rhs * scratch.phi_field[i]
+ time_step *
gamma
* scratch.phi_field[i])
*
scratch.finite_element_values.JxW(q);
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- if (cell->get_fe().system_to_component_index(j).first == solution_component)
+ for (unsigned int j=0; j<advection_dofs_per_cell; ++j)
{
data.local_matrix(i,j)
+= (
@@ -1344,6 +1397,8 @@
Simulator<dim>::
copy_local_to_global_advection_system (const internal::Assembly::CopyData::AdvectionSystem<dim> &data)
{
+ // copy entries into the global matrix. note that these local contributions
+ // only correspond to the advection dofs, as assembled above
current_constraints.distribute_local_to_global (data.local_matrix,
data.local_rhs,
data.local_dof_indices,
@@ -1410,7 +1465,13 @@
// rounded up. do so. (note that x/2 rounded up
// equals (x+1)/2 using integer division.)
internal::Assembly::Scratch::
- AdvectionSystem<dim> (finite_element, mapping,
+ AdvectionSystem<dim> (finite_element,
+ finite_element.base_element(temperature_or_composition.is_temperature()
+ ?
+ introspection.block_indices.temperature
+ :
+ introspection.block_indices.compositional_fields[temperature_or_composition.compositional_variable]),
+ mapping,
QGauss<dim>((temperature_or_composition.is_temperature()
?
parameters.temperature_degree
@@ -1420,7 +1481,11 @@
(parameters.stokes_velocity_degree+1)/2),
parameters.n_compositional_fields),
internal::Assembly::CopyData::
- AdvectionSystem<dim> (finite_element));
+ AdvectionSystem<dim> (finite_element.base_element(temperature_or_composition.is_temperature()
+ ?
+ introspection.block_indices.temperature
+ :
+ introspection.block_indices.compositional_fields[temperature_or_composition.compositional_variable])));
system_matrix.compress(VectorOperation::add);
system_rhs.compress(VectorOperation::add);
More information about the CIG-COMMITS
mailing list