[cig-commits] commit 1998 by buerg to /var/svn/dealii/aspect
dealii.demon at gmail.com
dealii.demon at gmail.com
Wed Oct 30 09:13:48 PDT 2013
Revision 1998
Fix bug in solver.
U trunk/aspire/include/aspect/simulator.h
U trunk/aspire/source/simulator/assembly.cc
U trunk/aspire/source/simulator/core.cc
U trunk/aspire/source/simulator/solver.cc
http://www.dealii.org/websvn/revision.php?repname=Aspect+Repository&path=%2F&rev=1998&peg=1998
Diff:
Modified: trunk/aspire/include/aspect/simulator.h
===================================================================
--- trunk/aspire/include/aspect/simulator.h 2013-10-28 19:58:43 UTC (rev 1997)
+++ trunk/aspire/include/aspect/simulator.h 2013-10-30 16:12:59 UTC (rev 1998)
@@ -68,6 +68,7 @@
{
namespace Scratch
{
+ template <int dim> struct StokesPreconditioner;
template <int dim> struct StokesSystem;
template <int dim> struct TemperatureSystem;
template <int dim> struct CompositionSystem;
@@ -75,6 +76,7 @@
namespace CopyData
{
+ template <int dim> struct StokesPreconditioner;
template <int dim> struct StokesSystem;
template <int dim> struct TemperatureSystem;
template <int dim> struct CompositionSystem;
@@ -404,6 +406,14 @@
void advance_chemistry ();
/**
+ * Initiate the assembly of the Stokes preconditioner.
+ *
+ * This function is implemented in
+ * <code>source/simulator/assembly.cc</code>.
+ */
+ void assemble_stokes_preconditioner ();
+
+ /**
* Initiate the assembly of the Stokes matrix and right hand side.
*
* This function is implemented in
@@ -541,6 +551,15 @@
*/
/**
* Set up the size and structure of the matrix used to store the
+ * preconditioner for the linear system for the velocity.
+ *
+ * This function is implemented in
+ * <code>source/simulator/core.cc</code>.
+ */
+ void setup_preconditioner_velocity (const std::vector<IndexSet> &system_partitioning);
+
+ /**
+ * Set up the size and structure of the matrix used to store the
* elements of the linear system for the velocity.
*
* This function is implemented in
@@ -576,6 +595,19 @@
*/
/**
+ * Compute the integrals for the Stokes preconditioner on a single
+ * cell.
+ *
+ * This function is implemented in
+ * <code>source/simulator/assembly.cc</code>.
+ */
+ void
+ local_assemble_stokes_preconditioner (const SynchronousIterators<std_cxx1x::tuple<FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>,
+ FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> > > &cell,
+ internal::Assembly::Scratch::StokesPreconditioner<dim> &scratch,
+ internal::Assembly::CopyData::StokesPreconditioner<dim> &data);
+
+ /**
* Compute the integrals for the Stokes matrix and right hand side
* on a single cell.
*
@@ -589,6 +621,16 @@
internal::Assembly::CopyData::StokesSystem<dim> &data);
/**
+ * Copy the contribution to the Stokes preconditioner from a single
+ * cell into the preconditioner matrix that stores these elements.
+ *
+ * This function is implemented in
+ * <code>source/simulator/assembly.cc</code>.
+ */
+ void
+ copy_local_to_global_stokes_preconditioner (const internal::Assembly::CopyData::StokesPreconditioner<dim> &data);
+
+ /**
* Copy the contribution to the Stokes system
* from a single cell into the global matrix that stores these elements.
*
@@ -937,6 +979,7 @@
* @name Variables that describe the linear systems and solution vectors
* @{
*/
+ LinearAlgebra::BlockSparseMatrix preconditioner_velocity;
LinearAlgebra::BlockSparseMatrix system_matrix_velocity;
LinearAlgebra::SparseMatrix system_matrix_temperature;
LinearAlgebra::BlockSparseMatrix system_matrix_composition;
Modified: trunk/aspire/source/simulator/assembly.cc
===================================================================
--- trunk/aspire/source/simulator/assembly.cc 2013-10-28 19:58:43 UTC (rev 1997)
+++ trunk/aspire/source/simulator/assembly.cc 2013-10-30 16:12:59 UTC (rev 1998)
@@ -44,17 +44,17 @@
namespace Scratch
{
template <int dim>
- struct StokesSystem
+ struct StokesPreconditioner
{
- StokesSystem (const FiniteElement<dim> &fe_velocity,
- const FiniteElement<dim> &fe_advection,
- const Mapping<dim> &mapping,
- const Quadrature<dim> &quadrature,
- const UpdateFlags update_flags_velocity,
- const UpdateFlags update_flags_advection,
- const unsigned int n_compositional_fields);
+ StokesPreconditioner (const FiniteElement<dim> &fe_velocity,
+ const FiniteElement<dim> &fe_advection,
+ const Mapping<dim> &mapping,
+ const Quadrature<dim> &quadrature,
+ const UpdateFlags update_flags_velocity,
+ const UpdateFlags update_flags_advection,
+ const unsigned int n_compositional_fields);
- StokesSystem (const StokesSystem<dim> &data);
+ StokesPreconditioner (const StokesPreconditioner<dim> &data);
FEValues<dim> fe_values_velocity;
FEValues<dim> fe_values_advection;
@@ -64,13 +64,7 @@
std::vector<SymmetricTensor<2,dim> > symmetric_grads_phi_u;
std::vector<Tensor<2,dim> > grads_phi_u;
std::vector<double> div_phi_u;
- std::vector<double> temperature_values;
- std::vector<Tensor<1, dim> > temperature_gradients;
- std::vector<std::vector<double> > composition_values;
- std::vector<std::vector<Tensor<1, dim> > > composition_gradients;
- std::vector<Tensor<1, dim> > pressure_gradients;
std::vector<Tensor<1,dim> > velocity_values;
- std::vector<Tensor<1,dim> > old_velocity_values;
typename MaterialModel::Interface<dim>::MaterialModelInputs material_model_inputs;
typename MaterialModel::Interface<dim>::MaterialModelOutputs material_model_outputs;
@@ -79,14 +73,14 @@
template <int dim>
- StokesSystem<dim>::
- StokesSystem (const FiniteElement<dim> &fe_velocity,
- const FiniteElement<dim> &fe_advection,
- const Mapping<dim> &mapping,
- const Quadrature<dim> &quadrature,
- const UpdateFlags update_flags_velocity,
- const UpdateFlags update_flags_advection,
- const unsigned int n_compositional_fields)
+ StokesPreconditioner<dim>::
+ StokesPreconditioner (const FiniteElement<dim> &fe_velocity,
+ const FiniteElement<dim> &fe_advection,
+ const Mapping<dim> &mapping,
+ const Quadrature<dim> &quadrature,
+ const UpdateFlags update_flags_velocity,
+ const UpdateFlags update_flags_advection,
+ const unsigned int n_compositional_fields)
:
fe_values_velocity (mapping, fe_velocity, quadrature,
update_flags_velocity),
@@ -97,15 +91,7 @@
symmetric_grads_phi_u (fe_velocity.dofs_per_cell),
grads_phi_u (fe_velocity.dofs_per_cell),
div_phi_u (fe_velocity.dofs_per_cell),
- temperature_values (quadrature.size ()),
- temperature_gradients (quadrature.size ()),
- composition_values (n_compositional_fields,
- std::vector<double> (quadrature.size ())),
- composition_gradients (n_compositional_fields,
- std::vector<Tensor<1, dim> > (quadrature.size ())),
- pressure_gradients (quadrature.size ()),
velocity_values (quadrature.size()),
- old_velocity_values (quadrature.size()),
material_model_inputs (quadrature.size (), n_compositional_fields),
material_model_outputs (quadrature.size (), n_compositional_fields)
{}
@@ -113,8 +99,8 @@
template <int dim>
- StokesSystem<dim>::
- StokesSystem (const StokesSystem<dim> &scratch)
+ StokesPreconditioner<dim>::
+ StokesPreconditioner (const StokesPreconditioner<dim> &scratch)
:
fe_values_velocity (scratch.fe_values_velocity.get_mapping (),
scratch.fe_values_velocity.get_fe (),
@@ -129,18 +115,71 @@
symmetric_grads_phi_u (scratch.symmetric_grads_phi_u),
grads_phi_u (scratch.grads_phi_u),
div_phi_u (scratch.div_phi_u),
- temperature_values (scratch.temperature_values),
- temperature_gradients (scratch.temperature_gradients),
- composition_values (scratch.composition_values),
- composition_gradients (scratch.composition_gradients),
- pressure_gradients (scratch.pressure_gradients),
velocity_values (scratch.velocity_values),
- old_velocity_values (scratch.old_velocity_values),
material_model_inputs(scratch.material_model_inputs),
material_model_outputs(scratch.material_model_outputs)
{}
template <int dim>
+ struct StokesSystem: public StokesPreconditioner<dim>
+ {
+ StokesSystem (const FiniteElement<dim> &fe_velocity,
+ const FiniteElement<dim> &fe_advection,
+ const Mapping<dim> &mapping,
+ const Quadrature<dim> &quadrature,
+ const UpdateFlags update_flags_velocity,
+ const UpdateFlags update_flags_advection,
+ const unsigned int n_compositional_fields);
+
+ StokesSystem (const StokesSystem<dim> &data);
+
+ std::vector<Tensor<1, dim> > pressure_gradients;
+ std::vector<std::vector<double> > composition_values;
+ std::vector<std::vector<Tensor<1, dim> > > composition_gradients;
+ std::vector<double> temperature_values;
+ std::vector<Tensor<1, dim> > temperature_gradients;
+ };
+
+
+
+ template <int dim>
+ StokesSystem<dim>::
+ StokesSystem (const FiniteElement<dim> &fe_velocity,
+ const FiniteElement<dim> &fe_advection,
+ const Mapping<dim> &mapping,
+ const Quadrature<dim> &quadrature,
+ const UpdateFlags update_flags_velocity,
+ const UpdateFlags update_flags_advection,
+ const unsigned int n_compositional_fields)
+ :
+ StokesPreconditioner<dim> (fe_velocity, fe_advection, mapping,
+ quadrature, update_flags_velocity,
+ update_flags_advection,
+ n_compositional_fields),
+ pressure_gradients (quadrature.size ()),
+ composition_values (n_compositional_fields,
+ std::vector<double> (quadrature.size ())),
+ composition_gradients (n_compositional_fields,
+ std::vector<Tensor<1, dim> > (quadrature.size ())),
+ temperature_values (quadrature.size ()),
+ temperature_gradients (quadrature.size ())
+ {}
+
+
+
+ template <int dim>
+ StokesSystem<dim>::
+ StokesSystem (const StokesSystem<dim> &scratch)
+ :
+ StokesPreconditioner<dim> (scratch),
+ pressure_gradients (scratch.pressure_gradients),
+ composition_values (scratch.composition_values),
+ composition_gradients (scratch.composition_gradients),
+ temperature_values (scratch.temperature_values),
+ temperature_gradients (scratch.temperature_gradients)
+ {}
+
+ template <int dim>
struct TemperatureSystem
{
TemperatureSystem (const FiniteElement<dim> &fe_velocity,
@@ -709,7 +748,8 @@
{
computing_timer.enter_section (" Build Stokes preconditioner");
pcout << " Building Stokes preconditioner ..." << std::flush;
-
+ // first assemble the raw matrices necessary for the preconditioner
+ assemble_stokes_preconditioner ();
// then extract the other information necessary to build the
// AMG preconditioners for the A and M blocks
std::vector<std::vector<bool> > constant_modes;
@@ -725,17 +765,155 @@
TrilinosWrappers::PreconditionBlockwiseDirect::AdditionalData Amg_data;
Amg_data.overlap = 4;
- prec_Mp->initialize (system_matrix_velocity.block (1,1));
- prec_A->initialize (system_matrix_velocity.block (0,0), Amg_data);
+ prec_Mp->initialize (preconditioner_velocity.block (1,1));
+ prec_A->initialize (preconditioner_velocity.block (0,0), Amg_data);
pcout << std::endl;
computing_timer.exit_section();
}
+
template <int dim>
+ void Simulator<dim>::assemble_stokes_preconditioner ()
+ {
+ preconditioner_velocity = 0.0;
+
+ const QGauss<dim> quadrature_formula (parameters.stokes_velocity_degree + 1);
+
+ FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
+ cell_filter_velocity_begin (IteratorFilters::LocallyOwnedCell (), dof_handler_velocity.begin_active ());
+ FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
+ cell_filter_advection_begin (IteratorFilters::LocallyOwnedCell (), dof_handler_advection.begin_active ());
+ typedef std_cxx1x::tuple<FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>,
+ FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> >
+ iterators_tuple;
+ SynchronousIterators<iterators_tuple>
+ filtered_iterators_begin (iterators_tuple (cell_filter_velocity_begin,
+ cell_filter_advection_begin));
+ FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
+ cell_filter_velocity_end (IteratorFilters::LocallyOwnedCell (), dof_handler_velocity.end ());
+ FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>
+ cell_filter_advection_end (IteratorFilters::LocallyOwnedCell (), dof_handler_advection.end ());
+ SynchronousIterators<iterators_tuple>
+ filtered_iterators_end (iterators_tuple (cell_filter_velocity_end,
+ cell_filter_advection_end));
+
+ WorkStream::
+ run (filtered_iterators_begin,
+ filtered_iterators_end,
+ std_cxx1x::bind (&Simulator<dim>::
+ local_assemble_stokes_preconditioner,
+ this,
+ std_cxx1x::_1,
+ std_cxx1x::_2,
+ std_cxx1x::_3),
+ std_cxx1x::bind (&Simulator<dim>::
+ copy_local_to_global_stokes_preconditioner,
+ this,
+ std_cxx1x::_1),
+ internal::Assembly::Scratch::
+ StokesPreconditioner<dim> (fe_velocity, fe_advection,
+ mapping, quadrature_formula,
+ (update_values |
+ update_quadrature_points |
+ update_JxW_values | update_gradients),
+ update_values,
+ parameters.n_compositional_fields),
+ internal::Assembly::CopyData::
+ StokesPreconditioner<dim> (fe_velocity));
+
+ preconditioner_velocity.compress (VectorOperation::add);
+ }
+
+
+ template <int dim>
void
Simulator<dim>::
+ local_assemble_stokes_preconditioner (const SynchronousIterators<std_cxx1x::tuple<FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>,
+ FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> > > &cell,
+ internal::Assembly::Scratch::StokesPreconditioner<dim> &scratch,
+ internal::Assembly::CopyData::StokesPreconditioner<dim> &data)
+ {
+ const unsigned int dofs_per_cell = scratch.fe_values_velocity.get_fe ().dofs_per_cell;
+ const unsigned int n_q_points = scratch.fe_values_velocity.n_quadrature_points;
+
+ scratch.fe_values_velocity.reinit (std_cxx1x::get<0> (cell.iterators));
+ scratch.fe_values_advection.reinit (std_cxx1x::get<1> (cell.iterators));
+ data.local_matrix = 0.0;
+
+ compute_material_model_input_values (current_linearization_point_velocity,
+ current_linearization_point_temperature,
+ current_linearization_point_composition,
+ scratch.fe_values_velocity,
+ scratch.fe_values_advection,
+ scratch.material_model_inputs);
+
+ material_model->evaluate(scratch.material_model_inputs,scratch.material_model_outputs);
+
+ const bool axisymmetric_terms = parameters.use_axisymmetric_discretization;
+
+ scratch.fe_values_velocity[introspection.extractors.velocities].get_function_values (current_linearization_point_velocity,
+ scratch.velocity_values);
+
+ for (unsigned int q = 0; q < n_q_points; ++q)
+ {
+ for (unsigned int k = 0; k < dofs_per_cell; ++k)
+ {
+ scratch.phi_u[k] = scratch.fe_values_velocity[introspection.extractors.velocities].value (k, q);
+ scratch.phi_p[k] = scratch.fe_values_velocity[introspection.extractors.pressure].value (k, q);
+ scratch.grads_phi_u[k] = scratch.fe_values_velocity[introspection.extractors.velocities].gradient (k, q);
+ scratch.symmetric_grads_phi_u[k]
+ = scratch.fe_values_velocity[introspection.extractors.velocities].symmetric_gradient (k, q);
+ scratch.div_phi_u[k] = scratch.fe_values_velocity[introspection.extractors.velocities].divergence (k, q);
+ }
+
+ const double eta = scratch.material_model_outputs.viscosities[q];
+ const double density = scratch.material_model_outputs.densities[q];
+
+ 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) += ((2.0 * eta * (scratch.symmetric_grads_phi_u[i] * scratch.symmetric_grads_phi_u[j])
+ + density * (scratch.phi_u[i] * (scratch.grads_phi_u[j] * scratch.velocity_values[q]))
+ + (material_model->is_compressible()
+ ? -2.0/3.0 * eta + parameters.stabilization_graddiv
+ : parameters.stabilization_graddiv)
+ * (scratch.div_phi_u[i] * scratch.div_phi_u[j])
+ + scratch.phi_p[i] * scratch.phi_p[j] / (eta + parameters.stabilization_graddiv)
+ - 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)
+ + (axisymmetric_terms ?
+ 2.0 * eta * scratch.phi_u[i][0] * scratch.phi_u[j][0]
+ / scratch.fe_values_velocity.quadrature_point (q) (0)
+ - scratch.phi_u[i][0] * scratch.phi_p[j]
+ - scratch.phi_p[i] * scratch.phi_u[j][0] : 0))
+ * scratch.fe_values_velocity.JxW (q);
+ }
+ }
+
+ std_cxx1x::get<0> (cell.iterators)->get_dof_indices (data.local_dof_indices);
+ }
+
+
+
+ template <int dim>
+ void
+ Simulator<dim>::
+ copy_local_to_global_stokes_preconditioner (const internal::Assembly::CopyData::StokesPreconditioner<dim> &data)
+ {
+ current_constraints_velocity.distribute_local_to_global (data.local_matrix,
+ data.local_dof_indices,
+ preconditioner_velocity);
+ }
+
+
+ template <int dim>
+ void
+ Simulator<dim>::
local_assemble_stokes_system (const SynchronousIterators<std_cxx1x::tuple<FilteredIterator<typename DoFHandler<dim>::active_cell_iterator>,
FilteredIterator<typename DoFHandler<dim>::active_cell_iterator> > > &cell,
internal::Assembly::Scratch::StokesSystem<dim> &scratch,
@@ -751,9 +929,6 @@
if (material_model->is_compressible())
data.local_pressure_shape_function_integrals = 0;
-
- if (material_model->is_compressible())
- data.local_pressure_shape_function_integrals = 0;
compute_material_model_input_values (current_linearization_point_velocity,
current_linearization_point_temperature,
@@ -770,8 +945,6 @@
scratch.pressure_gradients);
scratch.fe_values_velocity[introspection.extractors.velocities].get_function_values (current_linearization_point_velocity,
scratch.velocity_values);
- scratch.fe_values_velocity[introspection.extractors.velocities].get_function_values (old_solution_velocity,
- scratch.old_velocity_values);
scratch.fe_values_advection.get_function_gradients (current_linearization_point_temperature,
scratch.temperature_gradients);
scratch.fe_values_advection.get_function_values (current_linearization_point_temperature,
@@ -791,11 +964,10 @@
{
scratch.phi_u[k] = scratch.fe_values_velocity[introspection.extractors.velocities].value (k, q);
scratch.phi_p[k] = scratch.fe_values_velocity[introspection.extractors.pressure].value (k, q);
- scratch.grads_phi_u[k]
- = scratch.fe_values_velocity[introspection.extractors.velocities].gradient(k,q);
+ scratch.grads_phi_u[k] = scratch.fe_values_velocity[introspection.extractors.velocities].gradient(k,q);
scratch.symmetric_grads_phi_u[k]
= scratch.fe_values_velocity[introspection.extractors.velocities].symmetric_gradient(k,q);
- scratch.div_phi_u[k] = scratch.fe_values_velocity[introspection.extractors.velocities].divergence (k, q);
+ scratch.div_phi_u[k] = scratch.fe_values_velocity[introspection.extractors.velocities].divergence (k, q);
}
const double eta = scratch.material_model_outputs.viscosities[q];
@@ -817,7 +989,6 @@
? -2.0/3.0 * eta + parameters.stabilization_graddiv
: parameters.stabilization_graddiv)
* (scratch.div_phi_u[i] * scratch.div_phi_u[j])
- + scratch.phi_p[i] * scratch.phi_p[j] / (eta + parameters.stabilization_graddiv)
- scratch.div_phi_u[i] * scratch.phi_p[j]
- scratch.phi_p[i] * scratch.div_phi_u[j])
*
@@ -866,12 +1037,8 @@
void Simulator<dim>::assemble_stokes_system ()
{
computing_timer.enter_section (" Assemble Stokes system");
- system_matrix_velocity.block (introspection.block_indices.pressure,
- introspection.block_indices.pressure) = 0;
- system_matrix_velocity.block (introspection.block_indices.velocities,
- introspection.block_indices.velocities) = 0;
- system_rhs_velocity.block (introspection.block_indices.pressure) = 0;
- system_rhs_velocity.block (introspection.block_indices.velocities) = 0;
+ system_matrix_velocity = 0.0;
+ system_rhs_velocity = 0.0;
if (material_model->is_compressible())
pressure_shape_function_integrals = 0;
@@ -1002,7 +1169,7 @@
const double c_P = scratch.material_model_outputs.specific_heat[q];
const double lambda = scratch.material_model_outputs.thermal_conductivities[q];
const double old_value = scratch.old_temperature_values[q];
- const double rho = scratch.material_model_outputs.densities[q];
+ const double rho = 0.0*scratch.material_model_outputs.densities[q];
const Tensor<1, dim> current_u = scratch.current_velocity_values[q];
for (unsigned int i = 0; i < dofs_per_cell; ++i)
@@ -1070,7 +1237,7 @@
scratch.phi_field[i] = scratch.fe_values_advection.shape_value (i, q);
}
- const double rho = scratch.material_model_outputs.densities[q];
+ const double rho = 0.0*scratch.material_model_outputs.densities[q];
const Tensor<1, dim> current_u = scratch.current_velocity_values[q];
for (unsigned int c = 0; c < parameters.n_compositional_fields; ++c)
Modified: trunk/aspire/source/simulator/core.cc
===================================================================
--- trunk/aspire/source/simulator/core.cc 2013-10-28 19:58:43 UTC (rev 1997)
+++ trunk/aspire/source/simulator/core.cc 2013-10-30 16:12:59 UTC (rev 1998)
@@ -541,9 +541,9 @@
template <int dim>
void
Simulator<dim>::
- setup_system_matrix_velocity (const std::vector<IndexSet> &system_partitioning)
+ setup_preconditioner_velocity (const std::vector<IndexSet> &system_partitioning)
{
- system_matrix_velocity.clear ();
+ preconditioner_velocity.clear ();
TrilinosWrappers::BlockSparsityPattern sp (system_partitioning,
mpi_communicator);
@@ -560,6 +560,38 @@
this_mpi_process(mpi_communicator));
sp.compress();
+ preconditioner_velocity.reinit (sp);
+ }
+
+
+
+ template <int dim>
+ void
+ Simulator<dim>::
+ setup_system_matrix_velocity (const std::vector<IndexSet> &system_partitioning)
+ {
+ system_matrix_velocity.clear ();
+
+ TrilinosWrappers::BlockSparsityPattern sp (system_partitioning,
+ mpi_communicator);
+ Table<2, DoFTools::Coupling> coupling (dim + 1, dim + 1);
+
+ for (unsigned int i = 0; i < dim; ++i)
+ {
+ for (unsigned int j = 0; j < dim; ++j)
+ coupling[i][j] = DoFTools::always;
+
+ coupling[i][dim] = DoFTools::always;
+ coupling[dim][i] = DoFTools::always;
+ }
+
+ DoFTools::make_sparsity_pattern (dof_handler_velocity,
+ coupling, sp,
+ constraints_velocity, false,
+ Utilities::MPI::
+ this_mpi_process(mpi_communicator));
+ sp.compress();
+
system_matrix_velocity.reinit (sp);
}
@@ -731,6 +763,7 @@
// finally initialize vectors, matrices, etc.
+ setup_preconditioner_velocity (introspection.index_sets.velocity_partitioning);
setup_system_matrix_velocity (introspection.index_sets.velocity_partitioning);
setup_system_matrix_temperature (introspection.index_sets.advection_partitioning);
setup_system_matrix_composition (introspection.index_sets.advection_partitioning);
Modified: trunk/aspire/source/simulator/solver.cc
===================================================================
--- trunk/aspire/source/simulator/solver.cc 2013-10-28 19:58:43 UTC (rev 1997)
+++ trunk/aspire/source/simulator/solver.cc 2013-10-30 16:12:59 UTC (rev 1998)
@@ -111,7 +111,7 @@
SolverControl solver_control(5000, 1e-6 * src.block(1).l2_norm());
TrilinosWrappers::SolverCG solver(solver_control);
-
+
// Trilinos reports a breakdown
// in case src=dst=0, even
// though it should return
@@ -252,7 +252,7 @@
distributed_stokes_solution.block (1) = current_linearization_point_velocity.block (1);
const double solver_tolerance = std::max (parameters.linear_stokes_solver_tolerance *
- system_rhs_velocity.l2_norm (), 1e-12);
+ system_rhs_velocity.l2_norm (), 1e-12) / parameters.damping;
PrimitiveVectorMemory<LinearAlgebra::BlockVector> mem;
SolverControl solver_control_cheap (parameters.n_cheap_stokes_solver_steps, solver_tolerance);
SolverControl solver_control_expensive (system_matrix_velocity.block (0, 1). m ()
@@ -269,7 +269,7 @@
// of only a single V-cycle
const internal::BlockSchurPreconditioner<TrilinosWrappers::PreconditionBlockwiseDirect,
LinearAlgebra::PreconditionILU>
- preconditioner (system_matrix_velocity, *Mp_preconditioner, *Amg_preconditioner, false);
+ preconditioner (preconditioner_velocity, *Mp_preconditioner, *Amg_preconditioner, false);
SolverFGMRES<LinearAlgebra::BlockVector> solver (solver_control_cheap, mem,
SolverFGMRES<LinearAlgebra::BlockVector>::AdditionalData (30, true));
@@ -282,7 +282,7 @@
{
const internal::BlockSchurPreconditioner<TrilinosWrappers::PreconditionBlockwiseDirect,
LinearAlgebra::PreconditionILU>
- preconditioner (system_matrix_velocity, *Mp_preconditioner, *Amg_preconditioner, true);
+ preconditioner (preconditioner_velocity, *Mp_preconditioner, *Amg_preconditioner, true);
SolverFGMRES<LinearAlgebra::BlockVector> solver (solver_control_expensive, mem,
SolverFGMRES<LinearAlgebra::BlockVector>::AdditionalData (50, true));
More information about the CIG-COMMITS
mailing list