Mfgmg

CMakeLists.txt

@@ -1,6 +1,7 @@
 # Set the name of the project and target:
 SET(TARGET "SAND")
 
+
 # Declare all source files the target consists of. Here, this is only
 # the one step-X.cc file, but as you expand your project you may wish
 # to add other source files as well. If your project becomes much larger,
@@ -21,9 +22,8 @@ SET(  TARGET_SRC ${TARGET_SRC} ${TARGET_INC}
 
 # Usually, you will not need to modify anything beyond this point...
 
-CMAKE_MINIMUM_REQUIRED(VERSION 2.8.12)
-
-FIND_PACKAGE(deal.II 9.2.0 QUIET
+CMAKE_MINIMUM_REQUIRED(VERSION 3.1.0)
+FIND_PACKAGE(deal.II 9.4.0 QUIET
   HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
   )
 IF(NOT ${deal.II_FOUND})

include/density_filter.h

@@ -20,10 +20,13 @@
 #include <deal.II/lac/packaged_operation.h>
 #include <deal.II/lac/sparse_direct.h>
 #include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/trilinos_parallel_block_vector.h>
+#include <deal.II/lac/generic_linear_algebra.h>
 
 #include <deal.II/grid/tria.h>
 #include <deal.II/grid/grid_generator.h>
 #include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/cell_id.h>
 
 #include <deal.II/dofs/dof_handler.h>
 #include <deal.II/dofs/dof_renumbering.h>
@@ -47,20 +50,37 @@
  * Once formed, we have F*\sigma = \rho*/
 namespace SAND {
     using namespace dealii;
-
+    namespace LA
+    {
+        using namespace dealii::LinearAlgebraTrilinos;
+    }
     template<int dim>
     class DensityFilter {
     public:
 
-        DensityFilter()=default;
+        MPI_Comm  mpi_communicator;
+        std::vector<IndexSet> owned_partitioning;
+        std::vector<IndexSet> relevant_partitioning;
 
-        SparseMatrix<double> filter_matrix;
+        DensityFilter();
+        DynamicSparsityPattern filter_dsp;
+        LA::MPI::SparseMatrix filter_matrix;
+        LA::MPI::SparseMatrix filter_matrix_transpose;
         SparsityPattern filter_sparsity_pattern;
-        void initialize(Triangulation<dim> &triangulation);
-
+        void initialize(DoFHandler<dim> &dof_handler);
+        std::set<types::global_dof_index> find_relevant_neighbors(types::global_dof_index cell_index) const;
 
     private:
-        std::set<typename Triangulation<dim>::cell_iterator> find_relevant_neighbors(typename Triangulation<dim>::cell_iterator cell) const;
+        std::vector<double> cell_m;
+        std::vector<double> x_coord;
+        std::vector<double> y_coord;
+        std::vector<double> z_coord;
+        std::vector<double> cell_m_part;
+        std::vector<double> x_coord_part;
+        std::vector<double> y_coord_part;
+        std::vector<double> z_coord_part;
+
+        ConditionalOStream pcout;
 
     };
 }

include/input_information.h

@@ -16,19 +16,23 @@ namespace SAND {
 
         //geometry options
         constexpr unsigned int geometry_base = GeometryOptions::mbb;
-        constexpr unsigned int dim = 2;
-        constexpr unsigned int refinements =  3;
+
+        constexpr unsigned int dim = 3;
+        // constexpr unsigned int refinements = 3;
+        extern unsigned int refinements;
 
         //nonlinear algorithm options
         constexpr double initial_barrier_size = 25;
-        constexpr double min_barrier_size = .00000;
-        constexpr double fraction_to_boundary = .9;
-        constexpr unsigned int max_steps=75;
+        constexpr double min_barrier_size = 1e-12;
+
+        constexpr double fraction_to_boundary = .7;
+        constexpr unsigned int max_steps=100;
+
         constexpr unsigned int barrier_reduction=BarrierOptions::loqo;
         constexpr double required_norm = .0001;
 
         //density filter options
-        constexpr double filter_r = .251;
+        constexpr double filter_r = .25;
 
         //other options
         constexpr double density_penalty_exponent = 3;
@@ -39,11 +43,22 @@ namespace SAND {
         constexpr bool output_parts_of_matrix = false;
 
         //Linear solver options
-        constexpr unsigned int solver_choice = SolverOptions::inexact_K_with_exact_A_gmres;
+        constexpr unsigned int solver_choice = SolverOptions::inexact_K_with_inexact_A_gmres;
         constexpr bool use_eisenstat_walker = false;
         constexpr double default_gmres_tolerance = 1e-6;
 
+        extern unsigned int a_inv_iterations;
+        extern unsigned int k_inv_iterations;
+
+        // constexpr double a_inv_iterations = 5;
+        // constexpr double k_inv_iterations = 30;
+
+        constexpr double a_rel_tol = 1e-12;
+        constexpr double k_rel_tol = 1e-12;
 
+        //Material Options
+        constexpr double material_lambda = 1.;
+        constexpr double material_mu = 1.;
     }
 }
-#endif //SAND_INPUT_INFORMATION_H
+#endif //SAND_INPUT_INFORMATION_H

include/kkt_system.h

@@ -5,11 +5,15 @@
 #ifndef SAND_KKT_SYSTEM_H
 #define SAND_KKT_SYSTEM_H
 #include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/conditional_ostream.h>
+#include <deal.II/base/index_set.h>
 #include <deal.II/base/function.h>
 #include <deal.II/base/tensor.h>
 #include <deal.II/base/timer.h>
 
 #include <deal.II/lac/block_vector.h>
+#include <deal.II/lac/vector.h>
 #include <deal.II/lac/full_matrix.h>
 #include <deal.II/lac/block_sparse_matrix.h>
 #include <deal.II/lac/solver_cg.h>
@@ -18,10 +22,18 @@
 #include <deal.II/lac/packaged_operation.h>
 #include <deal.II/lac/sparse_direct.h>
 #include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/generic_linear_algebra.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/sparsity_tools.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/lac/trilinos_block_sparse_matrix.h>
+#include <deal.II/lac/trilinos_parallel_block_vector.h>
+#include <deal.II/lac/trilinos_vector.h>
 
-#include <deal.II/grid/tria.h>
 #include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/manifold_lib.h>
 #include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/grid_tools.h>
 
 #include <deal.II/dofs/dof_handler.h>
 #include <deal.II/dofs/dof_renumbering.h>
@@ -38,19 +50,50 @@
 #include <deal.II/numerics/error_estimator.h>
 #include <deal.II/hp/fe_collection.h>
 
+#include <deal.II/distributed/tria.h>
+#include <deal.II/distributed/grid_refinement.h>
+
+#include <deal.II/multigrid/multigrid.h>
+#include <deal.II/multigrid/mg_transfer_matrix_free.h>
+#include <deal.II/multigrid/mg_tools.h>
+#include <deal.II/multigrid/mg_coarse.h>
+#include <deal.II/multigrid/mg_smoother.h>
+#include <deal.II/multigrid/mg_matrix.h>
+
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/operators.h>
+#include <deal.II/matrix_free/fe_evaluation.h>
+
+
 #include "../include/schur_preconditioner.h"
 #include "../include/density_filter.h"
+#include "matrix_free_elasticity.h"
+
+#include <deal.II/base/conditional_ostream.h>
 
 #include <iostream>
 #include <fstream>
 #include <algorithm>
 namespace SAND {
+
+
+
+    namespace LA
+    {
+        using namespace dealii::LinearAlgebraTrilinos;
+    }
     using namespace dealii;
 
     template<int dim>
     class KktSystem {
 
+    using SystemMFMatrixType = MF_Elasticity_Operator<dim, 1, double>;
+    using LevelMFMatrixType = MF_Elasticity_Operator<dim, 1, double>;
     public:
+        MPI_Comm  mpi_communicator;
+        std::vector<IndexSet> owned_partitioning;
+        std::vector<IndexSet> relevant_partitioning;
+
         KktSystem();
 
         void
@@ -66,63 +109,101 @@ namespace SAND {
         setup_block_system();
 
         void
-        assemble_block_system(const BlockVector<double> &state, const double barrier_size);
+        assemble_block_system(const LA::MPI::BlockVector &state, const double barrier_size);
 
-        BlockVector<double>
-        solve(const BlockVector<double> &state, double barrier_size);
+        LA::MPI::BlockVector
+        solve(const LA::MPI::BlockVector &state);
 
-        BlockVector<double>
+        LA::MPI::BlockVector
         get_initial_state();
 
         double
-        calculate_objective_value(const BlockVector<double> &state) const;
+        calculate_objective_value(const LA::MPI::BlockVector &state) const;
 
         double
-        calculate_barrier_distance(const BlockVector<double> &state) const;
+        calculate_barrier_distance(const LA::MPI::BlockVector &state) const;
 
         double
-        calculate_feasibility(const BlockVector<double> &state, const double barrier_size) const;
+        calculate_feasibility(const LA::MPI::BlockVector &state, const double barrier_size) const;
 
         double
-        calculate_convergence(const BlockVector<double> &state) const;
+        calculate_convergence(const LA::MPI::BlockVector &state) const;
 
         void
-        output(const BlockVector<double> &state, const unsigned int j) const;
+        output(const LA::MPI::BlockVector &state, const unsigned int j) const;
 
         void
         calculate_initial_rhs_error();
 
         double
-        calculate_rhs_norm(const BlockVector<double> &state, const double barrier_size) const;
+        calculate_rhs_norm(const LA::MPI::BlockVector &state, const double barrier_size) const;
 
         void
-        output_stl(const BlockVector<double> &state);
+        output_stl(const LA::MPI::BlockVector &state);
 
     private:
 
-        BlockVector<double>
-        calculate_rhs(const BlockVector<double> &test_solution, const double barrier_size) const;
+        LA::MPI::BlockVector
+        calculate_rhs(const LA::MPI::BlockVector &test_solution, const double barrier_size) const;
 
+        BlockDynamicSparsityPattern dsp;
         BlockSparsityPattern sparsity_pattern;
-        BlockSparseMatrix<double> system_matrix;
-        BlockVector<double> linear_solution;
-        BlockVector<double> system_rhs;
-        Triangulation<dim> triangulation;
+        mutable LA::MPI::BlockSparseMatrix system_matrix;
+        mutable LA::MPI::BlockVector locally_relevant_solution;
+        mutable LA::MPI::BlockVector distributed_solution;
+        LA::MPI::BlockVector system_rhs;
+        parallel::distributed::Triangulation<dim> triangulation;
         DoFHandler<dim> dof_handler;
+        DoFHandler<dim> dof_handler_displacement;
+        DoFHandler<dim> dof_handler_density;
+
+        std::map<types::global_dof_index,types::global_dof_index> displacement_to_system_dof_index_map;
+        MGLevelObject<std::map<types::global_dof_index,types::global_dof_index>> level_displacement_to_system_dof_index_map;
+
         AffineConstraints<double> constraints;
+        AffineConstraints<double> displacement_constraints;
         FESystem<dim> fe_nine;
         FESystem<dim> fe_ten;
         hp::FECollection<dim> fe_collection;
+        FESystem<dim> fe_displacement;
+        FE_DGQ<dim> fe_density;
         const double density_ratio;
         const double density_penalty_exponent;
 
-        DensityFilter<dim> density_filter;
+        mutable DensityFilter<dim> density_filter;
 
         std::map<types::global_dof_index, double> boundary_values;
-
+        MGLevelObject<std::map<types::global_dof_index, double>> level_boundary_values;
+        ConditionalOStream pcout;
 
         double initial_rhs_error;
 
+        MGConstrainedDoFs mg_constrained_dofs;
+        SystemMFMatrixType elasticity_matrix_mf;
+        MGLevelObject<LevelMFMatrixType> mg_matrices;
+
+        OperatorCellData<dim, GMGNumberType> active_cell_data;
+        MGLevelObject<OperatorCellData<dim, GMGNumberType>> level_cell_data;
+        dealii::LinearAlgebra::distributed::Vector<double> active_density_vector;
+        dealii::LinearAlgebra::distributed::Vector<double> relevant_density_vector;
+        MGLevelObject<dealii::LinearAlgebra::distributed::Vector<double>> level_density_vector;
+
+        MGTransferMatrixFree<dim,GMGNumberType> transfer;
+        MGTransferMatrixFree<dim, double> mg_transfer;
+
+        using SmootherType = PreconditionChebyshev<LevelMFMatrixType, LinearAlgebra::distributed::Vector<double>>;
+        mg::SmootherRelaxation<SmootherType,LinearAlgebra::distributed::Vector<double>> mg_smoother;
+        MGLevelObject<typename SmootherType::AdditionalData> smoother_data;
+        MGCoarseGridApplySmoother<LinearAlgebra::distributed::Vector<double>> mg_coarse;
+        MGLevelObject<std::set<types::boundary_id>> level_dirichlet_boundary_dofs;
+        MGLevelObject<AffineConstraints<double>> mg_level_constraints;
+        MGLevelObject<MatrixFreeOperators::MGInterfaceOperator<LevelMFMatrixType>> mg_interface_matrices;
+
+        std::set<types::boundary_id> dirichlet_boundary;
+
+        LinearAlgebra::distributed::Vector<double> distributed_displacement_sol;
+        LinearAlgebra::distributed::Vector<double> distributed_displacement_rhs;
+
     };
 }