Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Adding some helpful comments to GMG code #3780

Merged
merged 2 commits into from Aug 14, 2020
Merged

Adding some helpful comments to GMG code #3780

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
907f764
some comments
Aug 12, 2020
4f4a7b7
fixes
Aug 13, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
14 changes: 6 additions & 8 deletions include/aspect/stokes_matrix_free.h
View file
Open in desktop
Expand Up @@ -305,7 +305,10 @@ namespace aspect
virtual void evaluate_material_model()=0;

/**
* Add correction to system RHS for non-zero boundary condition.
* Add correction to system RHS for non-zero boundary condition. For more information
* on exactly what this correction is and why it is computed, see the deal.II tutorial
* step 50 section "LaplaceProblem::assemble_rhs()":
* https://www.dealii.org/developer/doxygen/deal.II/step_50.html#LaplaceProblemassemble_rhs
*/
tcclevenger marked this conversation as resolved.
Show resolved Hide resolved
virtual void correct_stokes_rhs()=0;

Expand Down Expand Up @@ -433,7 +436,8 @@ namespace aspect
void evaluate_material_model() override;

/**
* Add correction to system RHS for non-zero boundary condition.
* Add correction to system RHS for non-zero boundary condition. See description in
* StokesMatrixFreeHandler::correct_stokes_rhs() for more information.
*/
tcclevenger marked this conversation as resolved.
Show resolved Hide resolved
void correct_stokes_rhs() override;

Expand All @@ -443,12 +447,6 @@ namespace aspect
*/
void build_preconditioner() override;

/**
* Get the workload imbalance of the distribution
* of the level hierarchy.
*/
double get_workload_imbalance();
tcclevenger marked this conversation as resolved.
Show resolved Hide resolved

/**
tcclevenger marked this conversation as resolved.
Show resolved Hide resolved
* Declare parameters. (No actual parameters at the moment).
*/
Expand Down
130 changes: 60 additions & 70 deletions source/simulator/stokes_matrix_free.cc
View file
Open in desktop
Expand Up @@ -1237,40 +1237,6 @@ namespace aspect
}


template <int dim, int velocity_degree>
double StokesMatrixFreeHandlerImplementation<dim, velocity_degree>::get_workload_imbalance ()
{
unsigned int n_proc = Utilities::MPI::n_mpi_processes(sim.triangulation.get_communicator());
unsigned int n_global_levels = sim.triangulation.n_global_levels();

unsigned long long int work_estimate = 0;
unsigned long long int total_cells_in_hierarchy = 0;

for (int lvl=n_global_levels-1; lvl>=0; --lvl)
{
unsigned long long int work_estimate_this_level;
unsigned long long int total_cells_on_lvl;
unsigned long long int n_owned_cells_on_lvl = 0;

for (const auto &cell: sim.triangulation.cell_iterators_on_level(lvl))
if (cell->is_locally_owned_on_level())
n_owned_cells_on_lvl += 1;

work_estimate_this_level = dealii::Utilities::MPI::max(n_owned_cells_on_lvl,sim.triangulation.get_communicator());

work_estimate += work_estimate_this_level;

total_cells_on_lvl = dealii::Utilities::MPI::sum(n_owned_cells_on_lvl,sim.triangulation.get_communicator());

total_cells_in_hierarchy += total_cells_on_lvl;
}
double ideal_work = static_cast<double>(total_cells_in_hierarchy) / static_cast<double>(n_proc);
double workload_imbalance_ratio = work_estimate / ideal_work;

return workload_imbalance_ratio;
}


template <int dim, int velocity_degree>
void StokesMatrixFreeHandlerImplementation<dim, velocity_degree>::evaluate_material_model ()
{
Expand All @@ -1295,6 +1261,10 @@ namespace aspect
MaterialModel::MaterialModelInputs<dim> in(fe_values.n_quadrature_points, sim.introspection.n_compositional_fields);
MaterialModel::MaterialModelOutputs<dim> out(fe_values.n_quadrature_points, sim.introspection.n_compositional_fields);

// This function call computes a cellwise projection of data defined at quadrature points to
// a vector defined by the projection DoFHandler. As an input, we must define a lambda which returns
// a viscosity value for each quadrature point of the given cell. The projection is then stored in
// the active level viscosity vector provided.
Utilities::project_cellwise<dim, dealii::LinearAlgebra::distributed::Vector<double>>(*(sim.mapping),
dof_handler_projection,
0,
Expand Down Expand Up @@ -1346,12 +1316,13 @@ namespace aspect
update_values);
std::vector<double> values_on_quad;

// Create active mesh viscosity table. For DGQ0, this is one value per cell,
// for DGQ1 this is n_q_points values per cell.
// Create active mesh viscosity table.
{
const unsigned int n_cells = stokes_matrix.get_matrix_free()->n_macro_cells();
const unsigned int n_q_points = quadrature_formula.size();

// One value per cell is required for DGQ0 projection and n_q_points
// values per cell for DGQ1.
if (dof_handler_projection.get_fe().degree == 0)
active_viscosity_table.reinit(TableIndices<2>(n_cells, 1));
else if (dof_handler_projection.get_fe().degree == 1)
Expand All @@ -1376,11 +1347,13 @@ namespace aspect
&dof_handler_projection);
DG_cell->get_active_or_mg_dof_indices(local_dof_indices);

// For DGQ0, we simply use the viscosity at the single
// support point of the element. For DGQ1, we must project
// back to quadrature point values.
if (dof_handler_projection.get_fe().degree == 0)
active_viscosity_table(cell, 0)[i] = active_viscosity_vector(local_dof_indices[0]);
else
{
// For DGQ1, project back to quadrature point values
fe_values_projection.reinit(DG_cell);
fe_values_projection.get_function_values(active_viscosity_vector,
local_dof_indices,
Expand All @@ -1398,6 +1371,7 @@ namespace aspect

const bool is_compressible = sim.material_model->is_compressible();

// Store viscosity tables and other data into the active level matrix-free objects.
stokes_matrix.fill_cell_data(active_viscosity_table,
sim.pressure_scaling,
is_compressible);
Expand All @@ -1410,11 +1384,13 @@ namespace aspect
sim.pressure_scaling);
}

// Project active viscosity vector to multilevel vectors
const unsigned int n_levels = sim.triangulation.n_global_levels();
level_viscosity_vector = 0.;
level_viscosity_vector.resize(0,n_levels-1);

// Project the active level viscosity vector to multilevel vector representations
// using MG transfer objects. This transfer is based on the same linear operator used to
// transfer data inside a v-cycle.
MGTransferMatrixFree<dim,double> transfer;
transfer.build(dof_handler_projection);
transfer.interpolate_to_mg(dof_handler_projection,
Expand All @@ -1425,11 +1401,12 @@ namespace aspect

for (unsigned int level=0; level<n_levels; ++level)
{
// Create multilevel viscosity tables. For DGQ0, this is one value per cell,
// for DGQ1 this is n_q_points values per cell.
// Create viscosity tables on each level.
const unsigned int n_cells = mg_matrices_A_block[level].get_matrix_free()->n_macro_cells();
const unsigned int n_q_points = quadrature_formula.size();

// One value per cell is required for DGQ0 projection and n_q_points
// values per cell for DGQ1.
if (dof_handler_projection.get_fe().degree == 0)
level_viscosity_tables[level].reinit(TableIndices<2>(n_cells, 1));
else
Expand All @@ -1452,11 +1429,13 @@ namespace aspect
&dof_handler_projection);
DG_cell->get_active_or_mg_dof_indices(local_dof_indices);

// For DGQ0, we simply use the viscosity at the single
// support point of the element. For DGQ1, we must project
// back to quadrature point values.
if (dof_handler_projection.get_fe().degree == 0)
level_viscosity_tables[level](cell, 0)[i] = level_viscosity_vector[level](local_dof_indices[0]);
else
{
// For DGQ1, project back to quadrature point vaues
fe_values_projection.reinit(DG_cell);
fe_values_projection.get_function_values(level_viscosity_vector[level],
local_dof_indices,
Expand All @@ -1471,6 +1450,7 @@ namespace aspect
}
}

// Store viscosity tables and other data into the multigrid level matrix-free objects.
mg_matrices_A_block[level].fill_cell_data (level_viscosity_tables[level],
is_compressible);
mg_matrices_Schur_complement[level].fill_cell_data (level_viscosity_tables[level],
Expand All @@ -1491,6 +1471,7 @@ namespace aspect
stokes_matrix.initialize_dof_vector(rhs_correction);
stokes_matrix.initialize_dof_vector(u0);

// The vector u0 is a zero vector, but with correct boundary values.
u0 = 0;
rhs_correction = 0;
sim.current_constraints.distribute(u0);
Expand All @@ -1504,10 +1485,15 @@ namespace aspect
const bool use_viscosity_at_quadrature_points
= (active_viscosity_table.size(1) == velocity.n_q_points);

// Much like the matrix-free apply_add() functions compute a matrix-vector
// product by looping over cells and applying local matrix operations,
// here we apply the negative of the stokes_matrix operator to u0.
for (unsigned int cell=0; cell<stokes_matrix.get_matrix_free()->n_macro_cells(); ++cell)
{
VectorizedArray<double> viscosity_x_2 = 2.0*active_viscosity_table(cell, 0);

// We must use read_dof_values_plain() as to not overwrite boundary information
// with the zero boundary used by the stokes_matrix operator.
velocity.reinit (cell);
velocity.read_dof_values_plain (u0.block(0));
velocity.evaluate (false,true,false);
Expand Down Expand Up @@ -1546,6 +1532,7 @@ namespace aspect
}
rhs_correction.compress(VectorOperation::add);

// Copy to the correct vector type and add the correction to the system rhs.
LinearAlgebra::BlockVector stokes_rhs_correction (sim.introspection.index_sets.stokes_partitioning, sim.mpi_communicator);
internal::ChangeVectorTypes::copy(stokes_rhs_correction,rhs_correction);

Expand All @@ -1562,11 +1549,13 @@ namespace aspect
double final_linear_residual = numbers::signaling_nan<double>();

// Below we define all the objects needed to build the GMG preconditioner:
using vector_t = dealii::LinearAlgebra::distributed::Vector<double>;
using VectorType = dealii::LinearAlgebra::distributed::Vector<double>;

// ABlock GMG Smoother: Chebyshev, degree 4
using ASmootherType = PreconditionChebyshev<ABlockMatrixType,vector_t>;
mg::SmootherRelaxation<ASmootherType, vector_t>
// ABlock GMG Smoother: Chebyshev, degree 4. Parameter values were chosen
// by trial and error. We use a more powerful version of the smoother on the
// coarsest level than on the other levels.
using ASmootherType = PreconditionChebyshev<ABlockMatrixType,VectorType>;
mg::SmootherRelaxation<ASmootherType, VectorType>
mg_smoother_A;
{
MGLevelObject<typename ASmootherType::AdditionalData> smoother_data_A;
Expand All @@ -1590,9 +1579,11 @@ namespace aspect
mg_smoother_A.initialize(mg_matrices_A_block, smoother_data_A);
}

// Schur complement matrix GMG Smoother: Chebyshev, degree 4
using MSmootherType = PreconditionChebyshev<SchurComplementMatrixType,vector_t>;
mg::SmootherRelaxation<MSmootherType, vector_t>
// Schur complement matrix GMG Smoother: Chebyshev, degree 4. Parameter values
// were chosen by trial and error. We use a more powerful version of the smoother
// on the coarsest level than on the other levels.
using MSmootherType = PreconditionChebyshev<SchurComplementMatrixType,VectorType>;
mg::SmootherRelaxation<MSmootherType, VectorType>
mg_smoother_Schur(4);
{
MGLevelObject<typename MSmootherType::AdditionalData> smoother_data_Schur;
Expand All @@ -1616,17 +1607,15 @@ namespace aspect
mg_smoother_Schur.initialize(mg_matrices_Schur_complement, smoother_data_Schur);
}

// Estimate the eigenvalues for the Chebyshev smoothers. If not running with
// deal.II 9.2.0, the eigenvalue estimate will be performed during the first
// application of the Chebyshev smoother during the solve.
// Estimate the eigenvalues for the Chebyshev smoothers.

//TODO: The setup for the smoother (as well as the entire GMG setup) should
// be moved to an assembly timing block instead of the Stokes solve
// timing block (as is currently the case).
for (unsigned int level = 0; level<sim.triangulation.n_global_levels(); ++level)
{
vector_t temp_velocity;
vector_t temp_pressure;
VectorType temp_velocity;
VectorType temp_pressure;
mg_matrices_A_block[level].initialize_dof_vector(temp_velocity);
mg_matrices_Schur_complement[level].initialize_dof_vector(temp_pressure);

Expand All @@ -1638,11 +1627,11 @@ namespace aspect
// Coarse Solver is just an application of the Chebyshev smoother setup
// in such a way to be a solver
//ABlock GMG
MGCoarseGridApplySmoother<vector_t> mg_coarse_A;
MGCoarseGridApplySmoother<VectorType> mg_coarse_A;
mg_coarse_A.initialize(mg_smoother_A);

//Schur complement matrix GMG
MGCoarseGridApplySmoother<vector_t> mg_coarse_Schur;
MGCoarseGridApplySmoother<VectorType> mg_coarse_Schur;
mg_coarse_Schur.initialize(mg_smoother_Schur);

// Interface matrices
Expand All @@ -1651,38 +1640,38 @@ namespace aspect
mg_interface_matrices_A.resize(0, sim.triangulation.n_global_levels()-1);
for (unsigned int level=0; level<sim.triangulation.n_global_levels(); ++level)
mg_interface_matrices_A[level].initialize(mg_matrices_A_block[level]);
mg::Matrix<vector_t > mg_interface_A(mg_interface_matrices_A);
mg::Matrix<VectorType > mg_interface_A(mg_interface_matrices_A);

// Schur complement matrix GMG
MGLevelObject<MatrixFreeOperators::MGInterfaceOperator<SchurComplementMatrixType> > mg_interface_matrices_Schur;
mg_interface_matrices_Schur.resize(0, sim.triangulation.n_global_levels()-1);
for (unsigned int level=0; level<sim.triangulation.n_global_levels(); ++level)
mg_interface_matrices_Schur[level].initialize(mg_matrices_Schur_complement[level]);
mg::Matrix<vector_t > mg_interface_Schur(mg_interface_matrices_Schur);
mg::Matrix<VectorType > mg_interface_Schur(mg_interface_matrices_Schur);

// MG Matrix
mg::Matrix<vector_t > mg_matrix_A(mg_matrices_A_block);
mg::Matrix<vector_t > mg_matrix_Schur(mg_matrices_Schur_complement);
mg::Matrix<VectorType > mg_matrix_A(mg_matrices_A_block);
mg::Matrix<VectorType > mg_matrix_Schur(mg_matrices_Schur_complement);

// MG object
// ABlock GMG
Multigrid<vector_t > mg_A(mg_matrix_A,
mg_coarse_A,
mg_transfer_A_block,
mg_smoother_A,
mg_smoother_A);
Multigrid<VectorType > mg_A(mg_matrix_A,
mg_coarse_A,
mg_transfer_A_block,
mg_smoother_A,
mg_smoother_A);
mg_A.set_edge_matrices(mg_interface_A, mg_interface_A);

// Schur complement matrix GMG
Multigrid<vector_t > mg_Schur(mg_matrix_Schur,
mg_coarse_Schur,
mg_transfer_Schur_complement,
mg_smoother_Schur,
mg_smoother_Schur);
Multigrid<VectorType > mg_Schur(mg_matrix_Schur,
mg_coarse_Schur,
mg_transfer_Schur_complement,
mg_smoother_Schur,
mg_smoother_Schur);
mg_Schur.set_edge_matrices(mg_interface_Schur, mg_interface_Schur);

// GMG Preconditioner for ABlock and Schur complement
using GMGPreconditioner = PreconditionMG<dim, vector_t, MGTransferMatrixFree<dim,double> >;
using GMGPreconditioner = PreconditionMG<dim, VectorType, MGTransferMatrixFree<dim,double> >;
GMGPreconditioner prec_A(dof_handler_v, mg_A, mg_transfer_A_block);
GMGPreconditioner prec_Schur(dof_handler_p, mg_Schur, mg_transfer_Schur_complement);

Expand Down Expand Up @@ -2120,6 +2109,7 @@ namespace aspect
}

// Setup the matrix-free operators

// Stokes matrix
{
typename MatrixFree<dim,double>::AdditionalData additional_data;
Expand Down