[cig-commits] [commit] master: start work on DirectSolver by putting u+p into one block (2dc3d7f)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Mon May 19 20:13:15 PDT 2014
Repository : https://github.com/geodynamics/aspect
On branch : master
Link : https://github.com/geodynamics/aspect/compare/76d4275352ef2cae5de9a073acd1c03a92c2670c...4f3d06fd1f3754419813db37ec9ef7f0f6f3cb15
>---------------------------------------------------------------
commit 2dc3d7fc5ed3a0bd72433b078a821d6e7b275fed
Author: Timo Heister <timo.heister at gmail.com>
Date: Mon May 12 15:59:12 2014 -0400
start work on DirectSolver by putting u+p into one block
Still incomplete, have to fix many other places that assume specific
blocks. I had to rework Introspection::BlockIndices for this.
>---------------------------------------------------------------
2dc3d7fc5ed3a0bd72433b078a821d6e7b275fed
include/aspect/introspection.h | 14 +++++++++-----
include/aspect/simulator.h | 2 ++
source/simulator/core.cc | 36 +++++++++++++++++++++++++++++-------
source/simulator/introspection.cc | 20 ++++++++------------
source/simulator/parameters.cc | 1 +
source/simulator/solver.cc | 26 ++++++++++++++++++++++++++
6 files changed, 75 insertions(+), 24 deletions(-)
diff --git a/include/aspect/introspection.h b/include/aspect/introspection.h
index aa76549..9be8359 100644
--- a/include/aspect/introspection.h
+++ b/include/aspect/introspection.h
@@ -58,8 +58,11 @@ namespace aspect
* @param n_compositional_fields The number of compositional fields that
* will be used in this simulation. This is used in initializing the
* fields of this class.
+ * @param split_vel_pressure Set to true if velocity and pressure should
+ * be in separate blocks.
*/
- Introspection (const unsigned int n_compositional_fields);
+ Introspection (const unsigned int n_compositional_fields,
+ const bool split_vel_pressure);
/**
* @name Things that are independent of the current mesh
@@ -126,11 +129,12 @@ namespace aspect
*/
struct BlockIndices
{
- BlockIndices (const unsigned int n_compositional_fields);
+ BlockIndices (const unsigned int n_compositional_fields,
+ const bool split_vel_pressure);
- static const unsigned int velocities = 0;
- static const unsigned int pressure = 1;
- static const unsigned int temperature = 2;
+ const unsigned int velocities;
+ const unsigned int pressure;
+ const unsigned int temperature;
const std::vector<unsigned int> compositional_fields;
};
/**
diff --git a/include/aspect/simulator.h b/include/aspect/simulator.h
index aaf7ff8..322fc64 100644
--- a/include/aspect/simulator.h
+++ b/include/aspect/simulator.h
@@ -188,6 +188,8 @@ namespace aspect
unsigned int n_cheap_stokes_solver_steps;
double temperature_solver_tolerance;
double composition_solver_tolerance;
+
+ bool direct_stokes_solver;
/**
* @}
*/
diff --git a/source/simulator/core.cc b/source/simulator/core.cc
index c4ee8c6..7283751 100644
--- a/source/simulator/core.cc
+++ b/source/simulator/core.cc
@@ -86,7 +86,7 @@ namespace aspect
ParameterHandler &prm)
:
parameters (prm),
- introspection (parameters.n_compositional_fields),
+ introspection (parameters.n_compositional_fields, !parameters.direct_stokes_solver),
mpi_communicator (Utilities::MPI::duplicate_communicator (mpi_communicator_)),
pcout (std::cout,
(Utilities::MPI::
@@ -896,20 +896,42 @@ namespace aspect
IndexSet system_index_set = dof_handler.locally_owned_dofs();
introspection.index_sets.system_partitioning.clear ();
- introspection.index_sets.system_partitioning.push_back(system_index_set.get_view(0,n_u));
- introspection.index_sets.system_partitioning.push_back(system_index_set.get_view(n_u,n_u+n_p));
+ if (parameters.direct_stokes_solver)
+ {
+ introspection.index_sets.system_partitioning.push_back(system_index_set.get_view(0,n_u+n_p));
+ }
+ else
+ {
+ introspection.index_sets.system_partitioning.push_back(system_index_set.get_view(0,n_u));
+ introspection.index_sets.system_partitioning.push_back(system_index_set.get_view(n_u,n_u+n_p));
+ }
introspection.index_sets.system_partitioning.push_back(system_index_set.get_view(n_u+n_p,n_u+n_p+n_T));
introspection.index_sets.stokes_partitioning.clear ();
- introspection.index_sets.stokes_partitioning.push_back(system_index_set.get_view(0,n_u));
- introspection.index_sets.stokes_partitioning.push_back(system_index_set.get_view(n_u,n_u+n_p));
+ if (parameters.direct_stokes_solver)
+ {
+ introspection.index_sets.stokes_partitioning.push_back(system_index_set.get_view(0,n_u+n_p));
+ }
+ else
+ {
+ introspection.index_sets.stokes_partitioning.push_back(system_index_set.get_view(0,n_u));
+ introspection.index_sets.stokes_partitioning.push_back(system_index_set.get_view(n_u,n_u+n_p));
+ }
DoFTools::extract_locally_relevant_dofs (dof_handler,
introspection.index_sets.system_relevant_set);
introspection.index_sets.system_relevant_partitioning.clear ();
- introspection.index_sets.system_relevant_partitioning
+ if (parameters.direct_stokes_solver)
+ {
+ introspection.index_sets.system_relevant_partitioning
+ .push_back(introspection.index_sets.system_relevant_set.get_view(n_u,n_u+n_p));
+ }
+ else
+ {
+ introspection.index_sets.system_relevant_partitioning
.push_back(introspection.index_sets.system_relevant_set.get_view(0,n_u));
- introspection.index_sets.system_relevant_partitioning
+ introspection.index_sets.system_relevant_partitioning
.push_back(introspection.index_sets.system_relevant_set.get_view(n_u,n_u+n_p));
+ }
introspection.index_sets.system_relevant_partitioning
.push_back(introspection.index_sets.system_relevant_set.get_view(n_u+n_p, n_u+n_p+n_T));
diff --git a/source/simulator/introspection.cc b/source/simulator/introspection.cc
index a41b6a9..fc5724b 100644
--- a/source/simulator/introspection.cc
+++ b/source/simulator/introspection.cc
@@ -41,15 +41,6 @@ namespace aspect
const unsigned int
Introspection<dim>::ComponentIndices::temperature;
- template <int dim>
- const unsigned int
- Introspection<dim>::BlockIndices::pressure;
-
- template <int dim>
- const unsigned int
- Introspection<dim>::BlockIndices::temperature;
-
-
namespace
{
template <int dim>
@@ -71,13 +62,14 @@ namespace aspect
template <int dim>
- Introspection<dim>::Introspection(const unsigned int n_compositional_fields)
+ Introspection<dim>::Introspection(const unsigned int n_compositional_fields,
+ const bool split_vel_pressure)
:
n_components (dim+2+n_compositional_fields),
n_blocks (3+n_compositional_fields),
extractors (n_compositional_fields),
component_indices (n_compositional_fields),
- block_indices (n_compositional_fields),
+ block_indices (n_compositional_fields, split_vel_pressure),
components_to_blocks (component_to_block_mapping<dim>(n_components)),
system_dofs_per_block (n_blocks)
{}
@@ -107,8 +99,12 @@ namespace aspect
template <int dim>
Introspection<dim>::BlockIndices::
- BlockIndices (const unsigned int n_compositional_fields)
+ BlockIndices (const unsigned int n_compositional_fields,
+ const bool split_vel_pressure)
:
+ velocities(0),
+ pressure (split_vel_pressure?1:0),
+ temperature (split_vel_pressure?2:1),
compositional_fields (half_open_sequence(3, 3+n_compositional_fields))
{}
diff --git a/source/simulator/parameters.cc b/source/simulator/parameters.cc
index 304dbaf..b12e1ea 100644
--- a/source/simulator/parameters.cc
+++ b/source/simulator/parameters.cc
@@ -547,6 +547,7 @@ namespace aspect
Simulator<dim>::Parameters::
parse_parameters (ParameterHandler &prm)
{
+ direct_stokes_solver = false;
// first, make sure that the ParameterHandler parser agrees
// with the code in main() about the meaning of the "Dimension"
// parameter
diff --git a/source/simulator/solver.cc b/source/simulator/solver.cc
index fa281af..d9c5c57 100644
--- a/source/simulator/solver.cc
+++ b/source/simulator/solver.cc
@@ -377,6 +377,32 @@ namespace aspect
{
computing_timer.enter_section (" Solve Stokes system");
+ if (parameters.direct_stokes_solver)
+ {
+ LinearAlgebra::BlockVector distributed_stokes_solution (introspection.index_sets.stokes_partitioning, mpi_communicator);
+
+ 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
+ distributed_stokes_solution.block(1) *= pressure_scaling;
+
+ // then copy back the solution from the temporary (non-ghosted) vector
+ // into the ghosted one with all solution components
+ solution.block(0) = distributed_stokes_solution.block(0);
+
+ remove_nullspace(solution, distributed_stokes_solution);
+
+ normalize_pressure(solution);
+
+ computing_timer.exit_section();
+
+ return 0;
+ }
+
pcout << " Solving Stokes system... " << std::flush;
const internal::StokesBlock stokes_block(system_matrix);
More information about the CIG-COMMITS
mailing list