diff --git a/include/base_interface.h b/include/base_interface.h index b8521421..2ce37175 100644 --- a/include/base_interface.h +++ b/include/base_interface.h @@ -194,6 +194,7 @@ class BaseInterface : public ParameterAcceptor, public SimulatorAccess &scratch, CopyData &data) const; +#ifdef DEAL_II_WITH_TRILINOS /** * Compute linear operators needed by the problem: - @p system_op * represents the system matrix associated to the Newton's @@ -226,6 +227,14 @@ class BaseInterface : public ParameterAcceptor, public SimulatorAccess &system_op, LinearOperator &prec_op, LinearOperator &prec_op_finer) const; +#endif //DEAL_II_WITH_TRILINOS + +#ifdef DEAL_II_WITH_PETSC + virtual void compute_system_operators(const std::vector >, + LinearOperator &system_op, + LinearOperator &prec_op, + LinearOperator &prec_op_finer) const; +#endif //DEAL_II_WITH_PETSC /** * Compute linear operators needed by the problem. When using diff --git a/include/interfaces/poisson_problem.h b/include/interfaces/poisson_problem.h index c6df0dc3..34889422 100644 --- a/include/interfaces/poisson_problem.h +++ b/include/interfaces/poisson_problem.h @@ -28,9 +28,14 @@ class PoissonProblem : public PDESystemInterface &, LinearOperator &) const; + void compute_system_operators(const std::vector >, + LinearOperator &, + LinearOperator &, + LinearOperator &) const; + private: mutable shared_ptr preconditioner; - + mutable shared_ptr PETSc_preconditioner; }; template @@ -120,4 +125,34 @@ PoissonProblem::compute_system_operators(const std::vector +void +PoissonProblem::compute_system_operators(const std::vector > matrices, + LinearOperator &system_op, + LinearOperator &prec_op, + LinearOperator &) const +{ + + PETSc_preconditioner.reset (new PETScWrappers::PreconditionJacobi()); + PETSc_preconditioner->initialize(matrices[0]->block(0,0)); + // + auto A = linear_operator( matrices[0]->block(0,0) ); + // + // LinearOperator P_inv; + // + auto Id = identity_operator(A.reinit_range_vector); //linear_operator(matrices[0]->block(0,0), *PETSc_preconditioner); + // + // auto P00 = P_inv; + // + // // ASSEMBLE THE PROBLEM: + system_op = block_operator<1, 1, LAPETSc::VectorType>({{ + {{ A }} + } + }); + + prec_op = block_operator<1, 1, LAPETSc::VectorType>({{ + {{ Id }} , + } + }); +} #endif diff --git a/include/lac/lac_initializer.h b/include/lac/lac_initializer.h index e2e22eb0..c86f2e01 100644 --- a/include/lac/lac_initializer.h +++ b/include/lac/lac_initializer.h @@ -22,6 +22,7 @@ class ScopedLACInitializer comm(comm) {}; +#ifdef DEAL_II_WITH_TRILINOS /** * Initialize a non ghosted TrilinosWrappers::MPI::BlockVector. */ @@ -30,7 +31,42 @@ class ScopedLACInitializer v.reinit(owned, comm, fast); }; + /** + * Initialize a non ghosted PETScWrappers::MPI::BlockVector. + */ + void operator() (TrilinosWrappers::BlockSparseMatrix &M, TrilinosWrappers::BlockSparsityPattern &sp, bool fast=false) + { + M.reinit(sp); + }; +#endif // DEAL_II_WITH_TRILINOS + +#ifdef DEAL_II_WITH_PETSC + /** + * Initialize a non ghosted PETScWrappers::MPI::BlockVector. + */ + void operator() (PETScWrappers::MPI::BlockVector &v, bool fast=false) + { + v.reinit(owned, comm); + }; + + /** + * Initialize a non ghosted PETScWrappers::MPI::BlockVector. + */ + void operator() (PETScWrappers::MPI::BlockSparseMatrix &M, dealii::BlockDynamicSparsityPattern &sp, bool fast=false) + { + M.reinit(owned, relevant, sp, comm); + }; +#endif // DEAL_II_WITH_PETSC + + /** + * Initialize a non ghosted PETScWrappers::MPI::BlockVector. + */ + void operator() (BlockSparseMatrix &M, dealii::BlockSparsityPattern &sp, bool fast=false) + { + M.reinit(sp); + }; +#ifdef DEAL_II_WITH_TRILINOS /** * Initialize a ghosted TrilinosWrappers::MPI::BlockVector. */ @@ -38,6 +74,17 @@ class ScopedLACInitializer { v.reinit(owned, relevant, comm, fast); }; +#endif // DEAL_II_WITH_TRILINOS + +#ifdef DEAL_II_WITH_PETSC + /** + * Initialize a ghosted PETScWrappers::MPI::BlockVector. + */ + void ghosted(PETScWrappers::MPI::BlockVector &v, bool fast=false) + { + v.reinit(owned, relevant, comm); + }; +#endif // DEAL_II_WITH_PETSC /** * Initialize a serial BlockVector. @@ -58,6 +105,7 @@ class ScopedLACInitializer (void)fast; }; +#ifdef DEAL_II_WITH_TRILINOS /** * Initialize a Trilinos Sparsity Pattern. */ @@ -74,6 +122,7 @@ class ScopedLACInitializer Utilities::MPI::this_mpi_process(comm)); s.compress(); } +#endif // DEAL_II_WITH_TRILINOS /** * Initialize a Deal.II Sparsity Pattern. @@ -93,6 +142,24 @@ class ScopedLACInitializer s.copy_from(csp); } + + /** + * Initialize a Deal.II Dynamic Sparsity Pattern. + */ + template + void operator() (dealii::BlockDynamicSparsityPattern &s, + const DoFHandler &dh, + const ConstraintMatrix &cm, + const Table<2,DoFTools::Coupling> &coupling) + { + s.reinit(dofs_per_block, dofs_per_block); + + DoFTools::make_sparsity_pattern (dh, + coupling, s, + cm, false); + s.compress(); + } + private: /** * Dofs per block. diff --git a/include/lac/lac_type.h b/include/lac/lac_type.h index 3715b3b4..505b8267 100644 --- a/include/lac/lac_type.h +++ b/include/lac/lac_type.h @@ -68,7 +68,7 @@ class LAPETSc */ typedef PETScWrappers::MPI::BlockSparseMatrix BlockMatrix; - typedef dealii::BlockSparsityPattern BlockSparsityPattern; + typedef dealii::BlockDynamicSparsityPattern BlockSparsityPattern; }; #endif // DEAL_II_WITH_PETSC diff --git a/include/pidomus.h b/include/pidomus.h index 9cbc9943..752e32c2 100644 --- a/include/pidomus.h +++ b/include/pidomus.h @@ -216,7 +216,7 @@ class piDoMUS : public ParameterAcceptor, public SundialsInterface &) const {} - +#ifdef DEAL_II_WITH_TRILINOS template void BaseInterface:: @@ -205,7 +205,20 @@ compute_system_operators(const std::vector > { Assert(false, ExcPureFunctionCalled ()); } +#endif //DEAL_II_WITH_TRILINOS +#ifdef DEAL_II_WITH_PETSC +template +void +BaseInterface:: +compute_system_operators(const std::vector >, + LinearOperator &, + LinearOperator &, + LinearOperator &) const +{ + Assert(false, ExcPureFunctionCalled ()); +} +#endif //DEAL_II_WITH_PETSC template void @@ -410,10 +423,18 @@ void BaseInterface::connect_to_signals() const {} +#ifdef DEAL_II_WITH_TRILINOS template class BaseInterface<2, 2, LATrilinos>; template class BaseInterface<2, 3, LATrilinos>; template class BaseInterface<3, 3, LATrilinos>; +#endif // DEAL_II_WITH_TRILINOS template class BaseInterface<2, 2, LADealII>; template class BaseInterface<2, 3, LADealII>; template class BaseInterface<3, 3, LADealII>; + +#ifdef DEAL_II_WITH_PETSC +template class BaseInterface<2, 2, LAPETSc>; +template class BaseInterface<2, 3, LAPETSc>; +template class BaseInterface<3, 3, LAPETSc>; +#endif // DEAL_II_WITH_PETSC \ No newline at end of file diff --git a/source/pidomus.cc b/source/pidomus.cc index d8715bc5..4bee7ddf 100644 --- a/source/pidomus.cc +++ b/source/pidomus.cc @@ -447,7 +447,7 @@ void piDoMUS::setup_dofs (const bool &first_run) *dof_handler, constraints, interface.get_matrix_coupling(i)); - matrices[i]->reinit(*matrix_sparsities[i]); + initializer(*matrices[i], *matrix_sparsities[i]); } if (first_run) @@ -681,6 +681,7 @@ void piDoMUS::assemble_matrices (const double t, /* ------------------------ MESH AND GRID ------------------------ */ +#ifdef DEAL_II_WITH_TRILINOS template void piDoMUS:: refine_and_transfer_solutions(LATrilinos::VectorType &y, @@ -742,6 +743,56 @@ refine_and_transfer_solutions(LATrilinos::VectorType &y, signals.end_refine_and_transfer_solutions(); } +#endif //DEAL_II_WITH_TRILINOS + +#ifdef DEAL_II_WITH_PETSC +template +void piDoMUS:: +refine_and_transfer_solutions(LAPETSc::VectorType &y, + LAPETSc::VectorType &y_dot, + LAPETSc::VectorType &y_expl, + LAPETSc::VectorType &distributed_y, + LAPETSc::VectorType &distributed_y_dot, + LAPETSc::VectorType &distributed_y_expl, + bool adaptive_refinement) +{ + distributed_y = y; + distributed_y_dot = y_dot; + distributed_y_expl = y_expl; + + parallel::distributed::SolutionTransfer > sol_tr(*dof_handler); + + std::vector old_sols (3); + old_sols[0] = &distributed_y; + old_sols[1] = &distributed_y_dot; + old_sols[2] = &distributed_y_expl; + + triangulation->prepare_coarsening_and_refinement(); + sol_tr.prepare_for_coarsening_and_refinement (old_sols); + + if (adaptive_refinement) + triangulation->execute_coarsening_and_refinement (); + else + triangulation->refine_global (1); + + setup_dofs(false); + + LAPETSc::VectorType new_sol (y); + LAPETSc::VectorType new_sol_dot (y_dot); + LAPETSc::VectorType new_sol_expl (y_expl); + + std::vector new_sols (3); + new_sols[0] = &new_sol; + new_sols[1] = &new_sol_dot; + new_sols[2] = &new_sol_expl; + + sol_tr.interpolate (new_sols); + + y = new_sol; + y_dot = new_sol_dot; + y_expl = new_sol_expl; +} +#endif //DEAL_II_WITH_PETSC template void piDoMUS:: @@ -1279,7 +1330,6 @@ piDoMUS::set_constrained_dofs_to_zero(typename LAC::VectorTy } } - template void piDoMUS::get_lumped_mass_matrix(typename LAC::VectorType &dst) const @@ -1351,11 +1401,18 @@ piDoMUS::get_lumped_mass_matrix(typename LAC::VectorType &ds } +#ifdef DEAL_II_WITH_TRILINOS template class piDoMUS<2, 2, LATrilinos>; template class piDoMUS<2, 3, LATrilinos>; template class piDoMUS<3, 3, LATrilinos>; +#endif // DEAL_II_WITH_TRILINOS template class piDoMUS<2, 2, LADealII>; template class piDoMUS<2, 3, LADealII>; template class piDoMUS<3, 3, LADealII>; +#ifdef DEAL_II_WITH_PETSC +template class piDoMUS<2, 2, LAPETSc>; +template class piDoMUS<2, 3, LAPETSc>; +template class piDoMUS<3, 3, LAPETSc>; +#endif // DEAL_II_WITH_PETSC diff --git a/source/simulator_access.cc b/source/simulator_access.cc index 1a373289..6b53b3f2 100644 --- a/source/simulator_access.cc +++ b/source/simulator_access.cc @@ -159,12 +159,18 @@ SimulatorAccess::get_fe () const } - +#ifdef DEAL_II_WITH_TRILINOS template class SimulatorAccess<2, 2, LATrilinos>; template class SimulatorAccess<2, 3, LATrilinos>; template class SimulatorAccess<3, 3, LATrilinos>; +#endif // DEAL_II_WITH_TRILINOS template class SimulatorAccess<2, 2, LADealII>; template class SimulatorAccess<2, 3, LADealII>; template class SimulatorAccess<3, 3, LADealII>; +#ifdef DEAL_II_WITH_PETSC +template class SimulatorAccess<2, 2, LAPETSc>; +template class SimulatorAccess<2, 3, LAPETSc>; +template class SimulatorAccess<3, 3, LAPETSc>; +#endif // DEAL_II_WITH_PETSC diff --git a/tests/parameters/poisson_problem_04.prm b/tests/parameters/poisson_problem_04.prm new file mode 100644 index 00000000..95d0f86a --- /dev/null +++ b/tests/parameters/poisson_problem_04.prm @@ -0,0 +1,150 @@ +# Parameter file generated with +# D2K_GIT_BRANCH= master +# D2K_GIT_SHORTREV= a43b2c9 +# DEAL_II_GIT_BRANCH= master +# DEAL_II_GIT_SHORTREV= 2aed525 +subsection Dirichlet boundary conditions + set IDs and component masks = 0=ALL + set IDs and expressions = 0=(1-y)*y*sin(2*pi*(x-t)) + set Known component names = u + set Used constants = +end +subsection Domain + set Colorize = false + set Copy boundary to manifold ids = false + set Copy material to manifold ids = false + set Create default manifolds = true + set Grid to generate = rectangle + set Input grid file name = + set Manifold descriptors = + set Mesh smoothing alogrithm = none + set Optional Point 1 = 0,0 + set Optional Point 2 = 1,1 + set Optional double 1 = 1.0 + set Optional double 2 = 0.5 + set Optional double 3 = 1.5 + set Optional int 1 = 1 + set Optional int 2 = 2 + set Optional vector of dim int = 1,1 + set Output grid file name = +end +subsection Error Tables + set Compute error = true + set Error file format = tex + set Error precision = 3 + set Output error tables = true + set Solution names = u + set Solution names for latex = u + set Table names = error + set Write error files = false + subsection Table 0 + set Add convergence rates = true + set Extra terms = cells,dofs + set Latex table caption = error + set List of error norms to compute = L2,H1 + set Rate key = + end +end +subsection Exact solution + set Function constants = + set Function expression = (1-y)*y*sin(2*pi*(x-t)) + set Variable names = x,y,t +end +subsection Forcing terms + set IDs and component masks = 0=u + set IDs and expressions = 0=2*pi*(y-1)*y*cos(-2*pi*(t-x))-2*(2*pi^2*(y-1)*y*sin(-2*pi*(t-x))-sin(-2*pi*(t-x))) + set Known component names = u + set Used constants = +end +subsection IDA Solver Parameters + set Absolute error tolerance = 1e-4 + set Final time = 0.1 + set Ignore algebraic terms for error computations = false + set Initial condition Newton max iterations = 5 + set Initial condition Newton parameter = 0.33 + set Initial condition type = use_y_diff + set Initial condition type after restart = use_y_dot + set Initial step size = 1e-4 + set Initial time = 0 + set Maximum number of nonlinear iterations = 10 + set Maximum order of BDF = 5 + set Min step size = 5e-5 + set Relative error tolerance = 1e-3 + set Seconds between each output = 1e-2 + set Show output of time steps = true + set Use local tolerances = false +end +subsection IMEX Parameters + set Absolute error tolerance = 1e-6 + set Final time = 0. + set Initial time = 0. + set Intervals between outputs = 1 + set Maximum number of inner nonlinear iterations = 3 + set Maximum number of outer nonlinear iterations = 5 + set Newton relaxation parameter = 1.000000 + set Relative error tolerance = 0.000000 + set Step size = 1e10 + set Update continuously Jacobian = true +end +subsection Initial solution + set Function constants = + set Function expression = 0 + set Variable names = x,y,t +end +subsection Initial solution_dot + set Function constants = + set Function expression = 0 + set Variable names = x,y,t +end +subsection Neumann boundary conditions + set IDs and component masks = + set IDs and expressions = + set Known component names = u + set Used constants = +end +subsection Output Parameters + set Files to save in run directory = + set Incremental run prefix = + set Output format = vtu + set Output partitioning = false + set Problem base name = strong + set Solution names = u + set Subdivisions = 1 +end +subsection Poisson problem + set Block of differential components = 1 + set Blocking of the finite element = u + set Finite element space = FESystem[FE_Q(1)] +end +subsection Refinement + set Bottom fraction = 0.2 + set Maximum number of cells (if available) = 0 + set Order (optimize) = 2 + set Refinement strategy = fraction + set Top fraction = 0.2 +end +subsection Time derivative of Dirichlet boundary conditions + set IDs and component masks = + set IDs and expressions = + set Known component names = u + set Used constants = +end +subsection Zero average constraints + set Known component names = u + set Zero average on boundary = + set Zero average on whole domain = +end +subsection pidomus + set Adaptive refinement = true + set Initial global refinement = 4 + set Jacobian solver tolerance = 1e-8 + set Maximum number of time steps = 10000 + set Number of cycles = 1 + set Overwrite Newton's iterations = true + set Print some useful informations about processes = true + set Refine mesh during transient = true + set Threshold for solver's restart = 1e-2 + set Time stepper = ida + set Show timer = false + set Use direct solver if available = true +end \ No newline at end of file diff --git a/tests/poisson_04.cc b/tests/poisson_04.cc new file mode 100644 index 00000000..fbb17374 --- /dev/null +++ b/tests/poisson_04.cc @@ -0,0 +1,29 @@ +#include "pidomus.h" +#include "interfaces/poisson_problem.h" +#include "tests.h" + +using namespace dealii; +int main (int argc, char *argv[]) +{ + + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + + initlog(); + deallog.depth_file(1); + + const int dim = 2; + const int spacedim = 2; + + PoissonProblem p; + piDoMUS solver ("pidomus",p); + ParameterAcceptor::initialize(SOURCE_DIR "/parameters/poisson_problem_04.prm", "used_parameters.prm"); + + + solver.run (); + + auto sol = solver.get_solution(); + for (unsigned int i = 0; i