From a7bdf6c8090cc7af98d5a50bf357a54d5c57b635 Mon Sep 17 00:00:00 2001 From: alberto sartori Date: Thu, 14 Jul 2016 11:46:45 +0200 Subject: [PATCH 1/2] get rid of sundials interface [skip ci] --- examples/heat_equation/CMakeLists.txt | 3 +- examples/heat_equation/include/heat_ida.h | 63 +++++++++------------ examples/heat_equation/source/heat_ida.cc | 67 +++++++++++++++++------ examples/heat_equation/source/main.cc | 4 -- 4 files changed, 78 insertions(+), 59 deletions(-) diff --git a/examples/heat_equation/CMakeLists.txt b/examples/heat_equation/CMakeLists.txt index f7315266..16bdff81 100644 --- a/examples/heat_equation/CMakeLists.txt +++ b/examples/heat_equation/CMakeLists.txt @@ -40,7 +40,8 @@ SET(_main source/main.cc) # Normally you shouldn't need to change anything below. ############################################################ # Declare all source files the target consists of: -file(GLOB _files source/*cc) +file(GLOB _files source/*cc + include/*.h) # Don't compile the main file into the library IF(NOT "${_main}" STREQUAL "") diff --git a/examples/heat_equation/include/heat_ida.h b/examples/heat_equation/include/heat_ida.h index d9507814..18a65285 100644 --- a/examples/heat_equation/include/heat_ida.h +++ b/examples/heat_equation/include/heat_ida.h @@ -33,7 +33,6 @@ #include -#include #include #include #include @@ -56,69 +55,57 @@ using namespace deal2lkit; typedef TrilinosWrappers::MPI::Vector VEC; template -class Heat : public SundialsInterface, public ParameterAcceptor +class Heat : public ParameterAcceptor { public: - Heat (const MPI_Comm &comm); + Heat (const MPI_Comm comm); virtual void declare_parameters (ParameterHandler &prm); void run (); - /********************************************************* - * Public interface from SundialsInterface - *********************************************************/ - virtual shared_ptr + shared_ptr create_new_vector() const; - /** Returns the number of degrees of freedom. Pure virtual function. */ - virtual unsigned int n_dofs() const; + /** Returns the number of degrees of freedom.*/ + unsigned int n_dofs() const; /** This function is called at the end of each iteration step for * the ode solver. Once again, the conversion between pointers and * other forms of vectors need to be done inside the inheriting * class. */ - virtual void output_step(const double t, - const VEC &solution, - const VEC &solution_dot, - const unsigned int step_number, - const double h); + void output_step(const double t, + const VEC &solution, + const VEC &solution_dot, + const unsigned int step_number); /** This function will check the behaviour of the solution. If it * is converged or if it is becoming unstable the time integrator * will be stopped. If the convergence is not achived the * calculation will be continued. If necessary, it can also reset * the time stepper. */ - virtual bool solver_should_restart(const double t, - const unsigned int step_number, - const double h, - VEC &solution, - VEC &solution_dot); + bool solver_should_restart(const double t, + VEC &solution, + VEC &solution_dot); /** For dae problems, we need a residual function. */ - virtual int residual(const double t, - const VEC &src_yy, - const VEC &src_yp, - VEC &dst); + int residual(const double t, + const VEC &src_yy, + const VEC &src_yp, + VEC &dst); /** Setup Jacobian system and preconditioner. */ - virtual int setup_jacobian(const double t, - const VEC &src_yy, - const VEC &src_yp, - const VEC &residual, - const double alpha); + int setup_jacobian(const double t, + const VEC &src_yy, + const VEC &src_yp, + const double alpha); /** Inverse of the Jacobian vector product. */ - virtual int solve_jacobian_system(const double t, - const VEC &y, - const VEC &y_dot, - const VEC &residual, - const double alpha, - const VEC &src, - VEC &dst) const; + int solve_jacobian_system(const VEC &src, + VEC &dst) const; @@ -128,7 +115,7 @@ class Heat : public SundialsInterface, public ParameterAcceptor corresponding variable is a differential component, zero otherwise. */ - virtual VEC &differential_components() const; + VEC &differential_components() const; private: void refine_mesh (); @@ -144,7 +131,7 @@ class Heat : public SundialsInterface, public ParameterAcceptor void set_constrained_dofs_to_zero(VEC &v) const; - const MPI_Comm &comm; + const MPI_Comm comm; unsigned int initial_global_refinement; unsigned int max_time_iterations; @@ -193,7 +180,7 @@ class Heat : public SundialsInterface, public ParameterAcceptor ParsedSolver Ainv; - IDAInterface dae; + IDAInterface ida; IndexSet global_partitioning; IndexSet partitioning; diff --git a/examples/heat_equation/source/heat_ida.cc b/examples/heat_equation/source/heat_ida.cc index d3edfe78..79e546d1 100644 --- a/examples/heat_equation/source/heat_ida.cc +++ b/examples/heat_equation/source/heat_ida.cc @@ -42,10 +42,9 @@ #include template -Heat::Heat (const MPI_Comm &communicator) +Heat::Heat (const MPI_Comm communicator) : - SundialsInterface (communicator), - comm(communicator), + comm(Utilities::MPI::duplicate_communicator(communicator)), pcout (std::cout, (Utilities::MPI::this_mpi_process(comm) == 0)), @@ -76,7 +75,7 @@ Heat::Heat (const MPI_Comm &communicator) /* iter= */ 1000, /* reduction= */1e-8, linear_operator(jacobian_matrix) ), - dae(*this) + ida("IDA Solver Parameters", comm) {} template @@ -424,8 +423,7 @@ template void Heat::output_step(const double t, const VEC &solution, const VEC &solution_dot, - const unsigned int step_number, - const double /* h */ ) + const unsigned int step_number) { computing_timer.enter_section ("Postprocessing"); dirichlet_bcs.set_time(t); @@ -464,8 +462,6 @@ void Heat::output_step(const double t, template bool Heat::solver_should_restart (const double , - const unsigned int , - const double , VEC &solution, VEC &solution_dot) { @@ -554,7 +550,6 @@ template int Heat::setup_jacobian (const double t, const VEC &src_yy, const VEC &src_yp, - const VEC &, const double alpha) { computing_timer.enter_section (" Setup Jacobian"); @@ -570,12 +565,7 @@ int Heat::setup_jacobian (const double t, } template -int Heat::solve_jacobian_system (const double , - const VEC &, - const VEC &, - const VEC &, - const double , - const VEC &src, +int Heat::solve_jacobian_system (const VEC &src, VEC &dst) const { computing_timer.enter_section (" Solve system"); @@ -619,7 +609,52 @@ void Heat::run () constraints.distribute(solution); - dae.start_ode(solution, solution_dot, max_time_iterations); + ida.create_new_vector = [this]() ->shared_ptr + { + return this->create_new_vector(); + }; + ida.residual = [this](const double t, + const VEC &y, + const VEC &y_dot, + VEC &residual) ->int + { + return this->residual(t,y,y_dot,residual); + }; + + ida.setup_jacobian = [this](const double t, + const VEC &y, + const VEC &y_dot, + const double alpha) ->int + { + return this->setup_jacobian(t,y,y_dot,alpha); + }; + + ida.solver_should_restart = [this](const double t, + VEC &y, + VEC &y_dot) ->bool + { + return this->solver_should_restart(t,y,y_dot); + }; + + ida.solve_jacobian_system = [this](const VEC &rhs, + VEC &dst) ->int + { + return this->solve_jacobian_system(rhs,dst); + }; + + ida.output_step = [this](const double t, + const VEC &y, + const VEC &y_dot, + const unsigned int step_number) + { + this->output_step(t,y,y_dot,step_number); + }; + ida.differential_components = [this]() ->VEC & + { + return this->differential_components(); + }; + + ida.solve_dae(solution, solution_dot); eh.error_from_exact(*mapping, *dof_handler, distributed_solution, exact_solution); eh.output_table(pcout); diff --git a/examples/heat_equation/source/main.cc b/examples/heat_equation/source/main.cc index b0721d69..8485c70f 100644 --- a/examples/heat_equation/source/main.cc +++ b/examples/heat_equation/source/main.cc @@ -30,10 +30,6 @@ int main(int argc, char **argv) MPI_Comm comm = MPI_COMM_WORLD; - int numprocs = Utilities::MPI::n_mpi_processes(comm); - int myid = Utilities::MPI::this_mpi_process(comm); - - Heat<2> solver(comm); ParameterAcceptor::initialize("../source/heat_ida.prm", "used_parameters.prm"); From b35bc605e2a2ad84ba944858d8b1f13dce043006 Mon Sep 17 00:00:00 2001 From: alberto sartori Date: Thu, 14 Jul 2016 12:59:22 +0200 Subject: [PATCH 2/2] updated dynamic_stokes example [skip ci] --- examples/dynamic_stokes/CMakeLists.txt | 3 +- examples/dynamic_stokes/include/stokes_ida.h | 63 +++++-------- examples/dynamic_stokes/source/stokes_ida.cc | 98 +++++++++++++------- 3 files changed, 93 insertions(+), 71 deletions(-) diff --git a/examples/dynamic_stokes/CMakeLists.txt b/examples/dynamic_stokes/CMakeLists.txt index b6659f3e..8378a2eb 100644 --- a/examples/dynamic_stokes/CMakeLists.txt +++ b/examples/dynamic_stokes/CMakeLists.txt @@ -40,7 +40,8 @@ SET(_main source/main.cc) # Normally you shouldn't need to change anything below. ############################################################ # Declare all source files the target consists of: -file(GLOB _files source/*cc) +file(GLOB _files source/*cc + include/*.h) # Don't compile the main file into the library IF(NOT "${_main}" STREQUAL "") diff --git a/examples/dynamic_stokes/include/stokes_ida.h b/examples/dynamic_stokes/include/stokes_ida.h index 07d413c9..ca2e24e4 100644 --- a/examples/dynamic_stokes/include/stokes_ida.h +++ b/examples/dynamic_stokes/include/stokes_ida.h @@ -36,7 +36,6 @@ #include #include -#include #include #include #include @@ -60,69 +59,57 @@ using namespace deal2lkit; typedef TrilinosWrappers::MPI::BlockVector VEC; template -class Stokes : public SundialsInterface, public ParameterAcceptor +class Stokes : public ParameterAcceptor { public: - Stokes (const MPI_Comm &comm); + Stokes (const MPI_Comm comm); virtual void declare_parameters (ParameterHandler &prm); void run (); - /********************************************************* - * Public interface from SundialsInterface - *********************************************************/ - virtual shared_ptr + shared_ptr create_new_vector() const; - /** Returns the number of degrees of freedom. Pure virtual function. */ - virtual unsigned int n_dofs() const; + /** Returns the number of degrees of freedom.*/ + unsigned int n_dofs() const; /** This function is called at the end of each iteration step for * the ode solver. Once again, the conversion between pointers and * other forms of vectors need to be done inside the inheriting * class. */ - virtual void output_step(const double t, - const VEC &solution, - const VEC &solution_dot, - const unsigned int step_number, - const double h); + void output_step(const double t, + const VEC &solution, + const VEC &solution_dot, + const unsigned int step_number); /** This function will check the behaviour of the solution. If it * is converged or if it is becoming unstable the time integrator * will be stopped. If the convergence is not achived the * calculation will be continued. If necessary, it can also reset * the time stepper. */ - virtual bool solver_should_restart(const double t, - const unsigned int step_number, - const double h, - VEC &solution, - VEC &solution_dot); + bool solver_should_restart(const double t, + VEC &solution, + VEC &solution_dot); /** For dae problems, we need a residual function. */ - virtual int residual(const double t, - const VEC &src_yy, - const VEC &src_yp, - VEC &dst); + int residual(const double t, + const VEC &src_yy, + const VEC &src_yp, + VEC &dst); /** Setup Jacobian system and preconditioner. */ - virtual int setup_jacobian(const double t, - const VEC &src_yy, - const VEC &src_yp, - const VEC &residual, - const double alpha); + int setup_jacobian(const double t, + const VEC &src_yy, + const VEC &src_yp, + const double alpha); /** Inverse of the Jacobian vector product. */ - virtual int solve_jacobian_system(const double t, - const VEC &y, - const VEC &y_dot, - const VEC &residual, - const double alpha, - const VEC &src, - VEC &dst) const; + int solve_jacobian_system(const VEC &src, + VEC &dst) const; @@ -132,7 +119,7 @@ class Stokes : public SundialsInterface, public ParameterAcceptor corresponding variable is a differential component, zero otherwise. */ - virtual VEC &differential_components() const; + VEC &differential_components() const; private: void refine_mesh (); @@ -150,7 +137,7 @@ class Stokes : public SundialsInterface, public ParameterAcceptor void update_constraints(const double &t); - const MPI_Comm &comm; + const MPI_Comm comm; unsigned int initial_global_refinement; unsigned int max_time_iterations; @@ -205,7 +192,7 @@ class Stokes : public SundialsInterface, public ParameterAcceptor ParsedDataOut data_out; - IDAInterface dae; + IDAInterface ida; IndexSet global_partitioning; std::vector partitioning; diff --git a/examples/dynamic_stokes/source/stokes_ida.cc b/examples/dynamic_stokes/source/stokes_ida.cc index 9b245e65..0fa90956 100644 --- a/examples/dynamic_stokes/source/stokes_ida.cc +++ b/examples/dynamic_stokes/source/stokes_ida.cc @@ -40,6 +40,7 @@ #include #include #include +#include #include #include @@ -47,10 +48,9 @@ #include template -Stokes::Stokes (const MPI_Comm &communicator) +Stokes::Stokes (const MPI_Comm communicator) : - SundialsInterface (communicator), - comm(communicator), + comm(Utilities::MPI::duplicate_communicator(communicator)), pcout (std::cout, (Utilities::MPI::this_mpi_process(comm) == 0)), @@ -77,7 +77,7 @@ Stokes::Stokes (const MPI_Comm &communicator) dirichlet_dot("Dirichlet dot", 3, "u,u,p", "0=u"), data_out("Output Parameters", "vtu"), - dae(*this) + ida("IDA Solver Parameters", comm) {} template @@ -152,7 +152,7 @@ void Stokes::setup_dofs (const bool &first_run) dof_handler->distribute_dofs (*fe); DoFRenumbering::component_wise (*dof_handler, sub_blocks); - mapping = SP(new MappingQ(2)); + mapping = SP(new MappingQ(2)); dofs_per_block.clear(); dofs_per_block.resize(fe_builder.n_blocks()); @@ -200,7 +200,7 @@ void Stokes::setup_dofs (const bool &first_run) DoFTools::make_hanging_node_constraints (*dof_handler, constraints); - dirichlet_bcs.interpolate_boundary_values(*dof_handler, constraints); + dirichlet_bcs.interpolate_boundary_values(*mapping,*dof_handler, constraints); constraints.close (); constraints_dot.clear (); @@ -209,7 +209,7 @@ void Stokes::setup_dofs (const bool &first_run) DoFTools::make_hanging_node_constraints (*dof_handler, constraints_dot); - dirichlet_dot.interpolate_boundary_values(*dof_handler, constraints_dot); + dirichlet_dot.interpolate_boundary_values(*mapping,*dof_handler, constraints_dot); constraints_dot.close (); jacobian_matrix.clear(); @@ -219,10 +219,9 @@ void Stokes::setup_dofs (const bool &first_run) Table<2,DoFTools::Coupling> coupling (dim+1, dim+1); for (unsigned int c=0; c::setup_dofs (const bool &first_run) for (unsigned int d=0; d::setup_dofs (const bool &first_run) if (first_run) { - VectorTools::interpolate(*dof_handler, initial_solution, solution); - VectorTools::interpolate(*dof_handler, initial_solution_dot, solution_dot); + VectorTools::interpolate(*mapping,*dof_handler, initial_solution, solution); + VectorTools::interpolate(*mapping,*dof_handler, initial_solution_dot, solution_dot); } computing_timer.exit_section(); @@ -284,7 +286,7 @@ void Stokes::update_constraints (const double &t) DoFTools::make_hanging_node_constraints (*dof_handler, constraints); - dirichlet_bcs.interpolate_boundary_values(*dof_handler, constraints); + dirichlet_bcs.interpolate_boundary_values(*mapping,*dof_handler, constraints); constraints.close (); constraints_dot.clear(); @@ -292,7 +294,7 @@ void Stokes::update_constraints (const double &t) DoFTools::make_hanging_node_constraints (*dof_handler, constraints_dot); - dirichlet_dot.interpolate_boundary_values(*dof_handler, constraints_dot); + dirichlet_dot.interpolate_boundary_values(*mapping,*dof_handler, constraints_dot); constraints_dot.close (); } @@ -327,7 +329,6 @@ void Stokes::assemble_jacobian_matrix(const double t, FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); FullMatrix cell_prec (dofs_per_cell, dofs_per_cell); - Vector cell_rhs (dofs_per_cell); const FEValuesExtractors::Vector velocities (0); const FEValuesExtractors::Scalar pressure (dim); @@ -501,10 +502,6 @@ int Stokes::residual (const double t, const FEValuesExtractors::Vector u (0); const FEValuesExtractors::Scalar p (dim); - std::vector > phi_u(dofs_per_cell); - std::vector > grads_phi_u(dofs_per_cell); - std::vector div_phi_u(dofs_per_cell); - typename DoFHandler::active_cell_iterator cell = dof_handler->begin_active(), endc = dof_handler->end(); @@ -603,8 +600,7 @@ template void Stokes::output_step(const double t, const VEC &solution, const VEC &solution_dot, - const unsigned int step_number, - const double /* h */ ) + const unsigned int step_number) { computing_timer.enter_section ("Postprocessing"); @@ -638,8 +634,6 @@ void Stokes::output_step(const double t, template bool Stokes::solver_should_restart (const double t, - const unsigned int , - const double , VEC &solution, VEC &solution_dot) { @@ -737,7 +731,6 @@ template int Stokes::setup_jacobian (const double t, const VEC &src_yy, const VEC &src_yp, - const VEC &, const double alpha) { computing_timer.enter_section (" Setup Jacobian"); @@ -750,12 +743,7 @@ int Stokes::setup_jacobian (const double t, } template -int Stokes::solve_jacobian_system (const double , - const VEC &, - const VEC &, - const VEC &, - const double , - const VEC &src, +int Stokes::solve_jacobian_system (const VEC &src, VEC &dst) const { computing_timer.enter_section (" Solve system"); @@ -848,7 +836,53 @@ void Stokes::run () constraints.distribute(solution); constraints_dot.distribute(solution_dot); - dae.start_ode(solution, solution_dot, max_time_iterations); + ida.create_new_vector = [this]() ->shared_ptr + { + return this->create_new_vector(); + }; + ida.residual = [this](const double t, + const VEC &y, + const VEC &y_dot, + VEC &residual) ->int + { + return this->residual(t,y,y_dot,residual); + }; + + ida.setup_jacobian = [this](const double t, + const VEC &y, + const VEC &y_dot, + const double alpha) ->int + { + return this->setup_jacobian(t,y,y_dot,alpha); + }; + + ida.solver_should_restart = [this](const double t, + VEC &y, + VEC &y_dot) ->bool + { + return this->solver_should_restart(t,y,y_dot); + }; + + ida.solve_jacobian_system = [this](const VEC &rhs, + VEC &dst) ->int + { + return this->solve_jacobian_system(rhs,dst); + }; + + ida.output_step = [this](const double t, + const VEC &y, + const VEC &y_dot, + const unsigned int step_number) + { + this->output_step(t,y,y_dot,step_number); + }; + + ida.differential_components = [this]() ->VEC & + { + return this->differential_components(); + }; + + ida.solve_dae(solution, solution_dot); computing_timer.print_summary(); timer_outfile.close();