[cig-commits] [commit] master: pressure scaling with direct solver (0699e7a)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Mon May 19 20:14:16 PDT 2014
Repository : https://github.com/geodynamics/aspect
On branch : master
Link : https://github.com/geodynamics/aspect/compare/76d4275352ef2cae5de9a073acd1c03a92c2670c...4f3d06fd1f3754419813db37ec9ef7f0f6f3cb15
>---------------------------------------------------------------
commit 0699e7ad4e5d7b93889698c171131b42399bb615
Author: Timo Heister <timo.heister at gmail.com>
Date: Mon May 19 18:08:37 2014 -0400
pressure scaling with direct solver
>---------------------------------------------------------------
0699e7ad4e5d7b93889698c171131b42399bb615
include/aspect/introspection.h | 6 +++
source/simulator/core.cc | 94 +++++++++++++++++++++++++++++++++++++++---
source/simulator/solver.cc | 16 +++++--
3 files changed, 108 insertions(+), 8 deletions(-)
diff --git a/include/aspect/introspection.h b/include/aspect/introspection.h
index 5139ada..44f719d 100644
--- a/include/aspect/introspection.h
+++ b/include/aspect/introspection.h
@@ -249,6 +249,12 @@ namespace aspect
* two elements of system_partitioning.
*/
std::vector<IndexSet> stokes_partitioning;
+
+ /**
+ * Pressure unknowns that are locally owned. This IndexSet is needed
+ * if velocity and pressure end up in the same block.
+ */
+ IndexSet locally_owned_pressure;
};
/**
* A variable that contains index sets describing which of the globally
diff --git a/source/simulator/core.cc b/source/simulator/core.cc
index 09ede6d..2c845d4 100644
--- a/source/simulator/core.cc
+++ b/source/simulator/core.cc
@@ -272,10 +272,7 @@ namespace aspect
// to make velocities and pressures of roughly the same (numerical) size,
// and we may have to fix up the right hand side vector before solving for
// compressible models if there are no in-/outflow boundaries
- if (parameters.use_direct_stokes_solver)
- pressure_scaling = 1.0;
- else
- pressure_scaling = material_model->reference_viscosity() / geometry_model->length_scale();
+ pressure_scaling = material_model->reference_viscosity() / geometry_model->length_scale();
std::set<types::boundary_id> open_velocity_boundary_indicators
= geometry_model->get_used_boundary_indicators();
@@ -866,6 +863,88 @@ namespace aspect
}
+ template <int dim, int spacedim>
+ std::vector<unsigned char>
+ get_local_component_association (const FiniteElement<dim,spacedim> &fe,
+ const ComponentMask &component_mask)
+ {
+ std::vector<unsigned char> local_component_association (fe.dofs_per_cell,
+ (unsigned char)(-1));
+
+ // compute the component each local dof belongs to.
+ // if the shape function is primitive, then this
+ // is simple and we can just associate it with
+ // what system_to_component_index gives us
+ for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+ if (fe.is_primitive(i))
+ local_component_association[i] =
+ fe.system_to_component_index(i).first;
+ else
+ // if the shape function is not primitive, then either use the
+ // component of the first nonzero component corresponding
+ // to this shape function (if the component is not specified
+ // in the component_mask), or the first component of this block
+ // that is listed in the component_mask (if the block this
+ // component corresponds to is indeed specified in the component
+ // mask)
+ {
+ const unsigned int first_comp =
+ fe.get_nonzero_components(i).first_selected_component();
+
+ if ((fe.get_nonzero_components(i)
+ &
+ component_mask).n_selected_components(fe.n_components()) == 0)
+ local_component_association[i] = first_comp;
+ else
+ // pick the component selected. we know from the previous 'if'
+ // that within the components that are nonzero for this
+ // shape function there must be at least one for which the
+ // mask is true, so we will for sure run into the break()
+ // at one point
+ for (unsigned int c=first_comp; c<fe.n_components(); ++c)
+ if (component_mask[c] == true)
+ {
+ local_component_association[i] = c;
+ break;
+ }
+ }
+
+ Assert (std::find (local_component_association.begin(),
+ local_component_association.end(),
+ (unsigned char)(-1))
+ ==
+ local_component_association.end(),
+ ExcInternalError());
+
+ return local_component_association;
+ }
+
+ template <int dim>
+ IndexSet extract_component_subset(DoFHandler<dim> & dof_handler, const ComponentMask & component_mask)
+ {
+ std::vector<unsigned char> local_asoc =
+ get_local_component_association (dof_handler.get_fe(),
+ ComponentMask(dof_handler.get_fe().n_components(), true));
+
+ IndexSet ret(dof_handler.n_dofs());
+
+ unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
+ std::vector<types::global_dof_index> indices(dofs_per_cell);
+ for (typename DoFHandler<dim>::active_cell_iterator cell=dof_handler.begin_active();
+ cell!=dof_handler.end(); ++cell)
+ if (cell->is_locally_owned())
+ {
+ cell->get_dof_indices(indices);
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ if (component_mask[local_asoc[i]])
+ ret.add_index(indices[i]);
+ }
+
+ return ret;
+ }
+
+
+
template <int dim>
void Simulator<dim>::setup_introspection ()
{
@@ -913,7 +992,12 @@ namespace aspect
introspection.index_sets.system_partitioning.push_back(system_index_set.get_view(0,n_u));
if (n_p != 0)
{
- introspection.index_sets.system_partitioning.push_back(system_index_set.get_view(n_u,n_u+n_p));
+ introspection.index_sets.locally_owned_pressure = system_index_set.get_view(n_u,n_u+n_p);
+ introspection.index_sets.system_partitioning.push_back(introspection.index_sets.locally_owned_pressure);
+ }
+ else
+ {
+ introspection.index_sets.locally_owned_pressure = system_index_set & extract_component_subset(dof_handler, introspection.component_masks.pressure);
}
introspection.index_sets.system_partitioning.push_back(system_index_set.get_view(n_u+n_p,n_u+n_p+n_T));
diff --git a/source/simulator/solver.cc b/source/simulator/solver.cc
index 765fdeb..96f18e4 100644
--- a/source/simulator/solver.cc
+++ b/source/simulator/solver.cc
@@ -381,15 +381,25 @@ namespace aspect
{
LinearAlgebra::BlockVector distributed_stokes_solution (introspection.index_sets.stokes_partitioning, mpi_communicator);
+ Assert(introspection.block_indices.velocities == 0, ExcNotImplemented());
+
+ if (material_model->is_compressible ())
+ make_pressure_rhs_compatible(system_rhs);
+
SolverControl cn;
TrilinosWrappers::SolverDirect solver(cn);
solver.solve(system_matrix.block(0,0), distributed_stokes_solution.block(0), system_rhs.block(0));
current_constraints.distribute (distributed_stokes_solution);
- // now rescale the pressure back to real physical units
- // we use pressure_scaling = 1.0 with a direct solver
- // distributed_stokes_solution.block(1) *= pressure_scaling;
+ // now rescale the pressure back to real physical units:
+ for (unsigned int i=0;i< introspection.index_sets.locally_owned_pressure.n_elements(); ++i)
+ {
+ types::global_dof_index idx = introspection.index_sets.locally_owned_pressure.nth_index_in_set(i);
+
+ distributed_stokes_solution(idx) *= pressure_scaling;
+ }
+ distributed_stokes_solution.compress(VectorOperation::insert);
// then copy back the solution from the temporary (non-ghosted) vector
// into the ghosted one with all solution components
More information about the CIG-COMMITS
mailing list