[cig-commits] [commit] master: introduce parameter for direct solver, some cleanup (502bdd3)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Mon May 19 20:13:58 PDT 2014
Repository : https://github.com/geodynamics/aspect
On branch : master
Link : https://github.com/geodynamics/aspect/compare/76d4275352ef2cae5de9a073acd1c03a92c2670c...4f3d06fd1f3754419813db37ec9ef7f0f6f3cb15
>---------------------------------------------------------------
commit 502bdd3751d083ee1d10b4c64df8a113fbe82055
Author: Timo Heister <timo.heister at gmail.com>
Date: Fri May 16 13:48:17 2014 -0400
introduce parameter for direct solver, some cleanup
>---------------------------------------------------------------
502bdd3751d083ee1d10b4c64df8a113fbe82055
include/aspect/simulator.h | 2 +-
source/simulator/assembly.cc | 2 +-
source/simulator/core.cc | 35 +++++++++++------------------------
source/simulator/helper_functions.cc | 3 ++-
source/simulator/parameters.cc | 9 ++++++++-
source/simulator/solver.cc | 2 +-
6 files changed, 24 insertions(+), 29 deletions(-)
diff --git a/include/aspect/simulator.h b/include/aspect/simulator.h
index 6609a97..e84b20c 100644
--- a/include/aspect/simulator.h
+++ b/include/aspect/simulator.h
@@ -183,13 +183,13 @@ namespace aspect
double surface_pressure;
double adiabatic_surface_temperature;
unsigned int timing_output_frequency;
+ bool use_direct_stokes_solver;
double linear_stokes_solver_tolerance;
unsigned int max_nonlinear_iterations;
unsigned int n_cheap_stokes_solver_steps;
double temperature_solver_tolerance;
double composition_solver_tolerance;
- bool direct_stokes_solver;
/**
* @}
*/
diff --git a/source/simulator/assembly.cc b/source/simulator/assembly.cc
index 673e482..4db6b28 100644
--- a/source/simulator/assembly.cc
+++ b/source/simulator/assembly.cc
@@ -1007,7 +1007,7 @@ namespace aspect
if (rebuild_stokes_preconditioner == false)
return;
- if (parameters.direct_stokes_solver)
+ if (parameters.use_direct_stokes_solver)
return;
computing_timer.enter_section (" Build Stokes preconditioner");
diff --git a/source/simulator/core.cc b/source/simulator/core.cc
index 2c77e6b..7a83cf6 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, !parameters.direct_stokes_solver),
+ introspection (parameters.n_compositional_fields, !parameters.use_direct_stokes_solver),
mpi_communicator (Utilities::MPI::duplicate_communicator (mpi_communicator_)),
pcout (std::cout,
(Utilities::MPI::
@@ -272,7 +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.direct_stokes_solver)
+ if (parameters.use_direct_stokes_solver)
pressure_scaling = 1.0;
else
pressure_scaling = material_model->reference_viscosity() / geometry_model->length_scale();
@@ -896,12 +896,12 @@ namespace aspect
n_p = introspection.system_dofs_per_block[introspection.block_indices.pressure],
n_T = introspection.system_dofs_per_block[introspection.block_indices.temperature];
- Assert(!parameters.direct_stokes_solver ||
+ Assert(!parameters.use_direct_stokes_solver ||
(introspection.block_indices.velocities == introspection.block_indices.pressure),
ExcInternalError());
// only count pressure once if they are in the same block
- if (parameters.direct_stokes_solver)
+ if (introspection.block_indices.velocities == introspection.block_indices.pressure)
n_p = 0;
std::vector<types::global_dof_index> n_C (parameters.n_compositional_fields);
@@ -912,41 +912,28 @@ namespace aspect
IndexSet system_index_set = dof_handler.locally_owned_dofs();
introspection.index_sets.system_partitioning.clear ();
- if (parameters.direct_stokes_solver)
+ 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(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 ();
- 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));
+ if (n_p != 0)
{
- 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 ();
- if (parameters.direct_stokes_solver)
- {
- introspection.index_sets.system_relevant_partitioning
- .push_back(introspection.index_sets.system_relevant_set.get_view(0,n_u+n_p));
- }
- else
+ introspection.index_sets.system_relevant_partitioning
+ .push_back(introspection.index_sets.system_relevant_set.get_view(0,n_u));
+ if (n_p != 0)
{
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
.push_back(introspection.index_sets.system_relevant_set.get_view(n_u,n_u+n_p));
}
introspection.index_sets.system_relevant_partitioning
diff --git a/source/simulator/helper_functions.cc b/source/simulator/helper_functions.cc
index f550237..488e298 100644
--- a/source/simulator/helper_functions.cc
+++ b/source/simulator/helper_functions.cc
@@ -535,7 +535,8 @@ namespace aspect
if (parameters.pressure_normalization == "no")
return;
- // TODO: I need to think about how to implement this. Choose "no" for now.
+ // TODO: pressure normalization currently does not work if velocity and
+ // pressure are in the same block.
Assert(introspection.block_indices.velocities != introspection.block_indices.pressure,
ExcNotImplemented());
diff --git a/source/simulator/parameters.cc b/source/simulator/parameters.cc
index b12e1ea..1414c4f 100644
--- a/source/simulator/parameters.cc
+++ b/source/simulator/parameters.cc
@@ -210,6 +210,13 @@ namespace aspect
"The name of the directory into which all output files should be "
"placed. This may be an absolute or a relative path.");
+ prm.declare_entry ("Use Stokes direct solver", "false",
+ Patterns::Bool(),
+ "If set to true the linear system for the Stokes equation will "
+ "be solved using Trilinos klu, otherwise an iterative Schur "
+ "complement solver is used. The direct solver is only efficient "
+ "for small problems.");
+
prm.declare_entry ("Linear solver tolerance", "1e-7",
Patterns::Double(0,1),
"A relative tolerance up to which the linear Stokes systems in each "
@@ -547,7 +554,6 @@ 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
@@ -610,6 +616,7 @@ namespace aspect
adiabatic_surface_temperature = prm.get_double ("Adiabatic surface temperature");
pressure_normalization = prm.get("Pressure normalization");
+ use_direct_stokes_solver = prm.get_bool("Use Stokes direct solver");
linear_stokes_solver_tolerance= prm.get_double ("Linear solver tolerance");
n_cheap_stokes_solver_steps = prm.get_integer ("Number of cheap Stokes solver steps");
temperature_solver_tolerance = prm.get_double ("Temperature solver tolerance");
diff --git a/source/simulator/solver.cc b/source/simulator/solver.cc
index 8c087a9..765fdeb 100644
--- a/source/simulator/solver.cc
+++ b/source/simulator/solver.cc
@@ -377,7 +377,7 @@ namespace aspect
{
computing_timer.enter_section (" Solve Stokes system");
- if (parameters.direct_stokes_solver)
+ if (parameters.use_direct_stokes_solver)
{
LinearAlgebra::BlockVector distributed_stokes_solution (introspection.index_sets.stokes_partitioning, mpi_communicator);
More information about the CIG-COMMITS
mailing list