diff --git a/include/pidomus.h b/include/pidomus.h index 9cbc9943..f85e1933 100644 --- a/include/pidomus.h +++ b/include/pidomus.h @@ -50,6 +50,21 @@ #include "lac/lac_type.h" #include "lac/lac_initializer.h" + +#include + +#include +#include +#include + +#ifdef DEAL_II_WITH_ZLIB +# include +#endif + +#include +#include + + using namespace dealii; using namespace deal2lkit; using namespace pidomus; @@ -290,6 +305,10 @@ class piDoMUS : public ParameterAcceptor, public SundialsInterfacesource/checkpoint_restart.cc. + */ + void create_snapshot() const; + + void save_solutions_and_triangulation(const LADealII::VectorType &y, + const LADealII::VectorType &y_dot, + const LADealII::VectorType &y_expl, + const LADealII::VectorType &, + const LADealII::VectorType &, + const LADealII::VectorType &) const; + + void save_solutions_and_triangulation(const LATrilinos::VectorType &, + const LATrilinos::VectorType &, + const LATrilinos::VectorType &, + const LATrilinos::VectorType &locally_relevant_y, + const LATrilinos::VectorType &locally_relevant_y_dot, + const LATrilinos::VectorType &locally_relevant_y_expl) const; + + + /** + * Restore the state of this program from a set of files in the output + * directory. + * + * This function is implemented in + * source/checkpoint_restart.cc. + */ + void resume_from_snapshot(); + + void load_solutions(LADealII::VectorType &y, + LADealII::VectorType &y_expl, + LADealII::VectorType &y_dot); + + void load_solutions(LATrilinos::VectorType &y, + LATrilinos::VectorType &y_expl, + LATrilinos::VectorType &y_dot); + +public: + /** + * Save a number of variables using BOOST serialization mechanism. + * + * This function is implemented in + * source/checkpoint_restart.cc. + */ + template + void serialize (Archive &ar, const unsigned int /*version*/) + { + ar ¤t_time; + ar ¤t_alpha; + ar ¤t_dt; + ar &step_number; + ar ¤t_cycle; + std::cout << "bannnnnnnnnnnnnnnnnnnaaa"<; + + }; #endif diff --git a/source/pidomus.cc b/source/pidomus.cc index 295f1f5c..8546b660 100644 --- a/source/pidomus.cc +++ b/source/pidomus.cc @@ -57,6 +57,17 @@ #include #include +#include +#include + +#ifdef DEAL_II_WITH_ZLIB +# include +#endif + +#include +#include + + #include #include "lac/lac_initializer.h" @@ -89,6 +100,7 @@ piDoMUS::piDoMUS (const std::string &name, // IDA calls residual first with alpha = 0 current_alpha(0.0), current_dt(std::numeric_limits::quiet_NaN()), + step_number(0), previous_time(std::numeric_limits::quiet_NaN()), previous_dt(std::numeric_limits::quiet_NaN()), second_to_last_time(std::numeric_limits::quiet_NaN()), @@ -245,6 +257,18 @@ declare_parameters (ParameterHandler &prm) "false", Patterns::Bool()); + add_parameter( prm, + &resume_computation, + "Resume computation from a snapshot", + "false", + Patterns::Bool()); + + add_parameter( prm, + &snapshot_folder, + "Snapshot folder", + "snapshots", + Patterns::Anything()); + } @@ -861,21 +885,31 @@ template void piDoMUS::run () { interface.set_stepper(time_stepper); + bool load_snap = true; for (current_cycle = 0; current_cycle < n_cycles; ++current_cycle) { - if (current_cycle == 0) + if (resume_computation && load_snap) { - make_grid_fe(); - setup_dofs(true); + resume_from_snapshot(); + load_snap = false; } else - refine_mesh(); - + { + if (current_cycle == 0) + { + make_grid_fe(); + setup_dofs(true); + } + else + refine_mesh(); + } constraints.distribute(solution); constraints_dot.distribute(solution_dot); if (time_stepper == "ida") + { ida.start_ode(solution, solution_dot, max_time_iterations); + } else if (time_stepper == "euler") { current_alpha = euler.get_alpha(); @@ -924,8 +958,9 @@ piDoMUS::output_step(const double t, { auto _timer = computing_timer.scoped_timer ("Postprocessing"); + this->step_number = step_number; syncronize(t,solution,solution_dot); - + create_snapshot(); interface.output_solution(current_cycle, step_number); } @@ -1349,6 +1384,296 @@ piDoMUS::get_lumped_mass_matrix(typename LAC::VectorType &ds } +template +void piDoMUS::create_snapshot() const +{ + auto _timer = computing_timer.scoped_timer ("Create snapshot"); + unsigned int my_id = Utilities::MPI::this_mpi_process (comm); + + if (my_id == 0) + { + // if we have previously written a snapshot, then keep the last + // snapshot in case this one fails to save. Note: static variables + // will only be initialied once per model run. + static bool previous_snapshot_exists = (resume_computation == true && file_exists(snapshot_folder + "restart.mesh")); + + if (previous_snapshot_exists == true) + { + rename_file (snapshot_folder + "restart.mesh", + snapshot_folder + "restart.mesh.old"); + rename_file (snapshot_folder + "restart.mesh.info", + snapshot_folder + "restart.mesh.info.old"); + rename_file (snapshot_folder + "restart.resume.z", + snapshot_folder + "restart.resume.z.old"); + } + // from now on, we know that if we get into this + // function again that a snapshot has previously + // been written + previous_snapshot_exists = true; + } + +// save Triangulation and Solution vectors: + + save_solutions_and_triangulation(solution,explicit_solution,solution_dot, + locally_relevant_solution, + locally_relevant_explicit_solution, + locally_relevant_solution_dot); + +// save general information This calls the serialization functions on all +// processes but only writes to the restart file on process 0 + { + std::ostringstream oss; + + // serialize into a stringstream + boost::archive::binary_oarchive oa (oss); + oa << (*this); + + // compress with zlib and write to file on the root processor +#ifdef DEAL_II_WITH_ZLIB + if (my_id == 0) + { + uLongf compressed_data_length = compressBound (oss.str().length()); + std::vector compressed_data (compressed_data_length); + int err = compress2 ((Bytef *) &compressed_data[0], + &compressed_data_length, + (const Bytef *) oss.str().data(), + oss.str().length(), + Z_BEST_COMPRESSION); + (void)err; + Assert (err == Z_OK, ExcInternalError()); + + // build compression header + const uint32_t compression_header[4] + = { 1, /* number of blocks */ + (uint32_t)oss.str().length(), /* size of block */ + (uint32_t)oss.str().length(), /* size of last block */ + (uint32_t)compressed_data_length + }; /* list of compressed sizes of blocks */ + + std::ofstream f ((snapshot_folder + "restart.resume.z").c_str()); + f.write((const char *)compression_header, 4 * sizeof(compression_header[0])); + f.write((char *)&compressed_data[0], compressed_data_length); + } +#else + AssertThrow (false, + ExcMessage ("You need to have deal.II configured with the 'libz' " + "option to support checkpoint/restart, but deal.II " + "did not detect its presence when you called 'cmake'.")); +#endif + + } + pcout << "*** Snapshot created!" << std::endl << std::endl; +} + + +template +void +piDoMUS:: +save_solutions_and_triangulation(const LADealII::VectorType &y, + const LADealII::VectorType &y_dot, + const LADealII::VectorType &y_expl, + const LADealII::VectorType &, + const LADealII::VectorType &, + const LADealII::VectorType &) const +{ + std::vector x_system (3); + x_system[0] = y; + x_system[1] = y_dot; + x_system[2] = y_expl; + + + SolutionTransfer > + system_trans (*dof_handler); + + //system_trans.prepare_serialization (x_system); + + triangulation->save ((snapshot_folder + "restart.mesh").c_str()); +} + +template +void +piDoMUS:: +save_solutions_and_triangulation(const LATrilinos::VectorType &, + const LATrilinos::VectorType &, + const LATrilinos::VectorType &, + const LATrilinos::VectorType &locally_relevant_y, + const LATrilinos::VectorType &locally_relevant_y_dot, + const LATrilinos::VectorType &locally_relevant_y_expl) const +{ + std::vector x_system (3); + x_system[0] = &locally_relevant_y; + x_system[1] = &locally_relevant_y_expl; + x_system[2] = &locally_relevant_y_dot; + + + parallel::distributed::SolutionTransfer > + system_trans (*dof_handler); + + system_trans.prepare_serialization (x_system); + + triangulation->save ((snapshot_folder + "restart.mesh").c_str()); +} + + + +template +void piDoMUS::resume_from_snapshot() +{ + // first check existence of the two restart files + { + const std::string filename = snapshot_folder + "restart.mesh"; + std::ifstream in (filename.c_str()); + if (!in) + AssertThrow (false, + ExcMessage (std::string("You are trying to restart a previous computation, " + "but the restart file <") + + + filename + + + "> does not appear to exist!")); + } + { + const std::string filename = snapshot_folder + "restart.resume.z"; + std::ifstream in (filename.c_str()); + if (!in) + AssertThrow (false, + ExcMessage (std::string("You are trying to restart a previous computation, " + "but the restart file <") + + + filename + + + "> does not appear to exist!")); + } + + pcout << "*** Resuming from snapshot!" << std::endl << std::endl; + + try + { + triangulation = SP(pgg.distributed(comm)); + triangulation->load ((snapshot_folder + "restart.mesh").c_str()); + dof_handler = SP(new DoFHandler(*triangulation)); + fe = SP(interface.pfe()); + } + catch (...) + { + AssertThrow(false, ExcMessage("Cannot open snapshot mesh file or read the triangulation stored there.")); + } + setup_dofs(false); + + load_solutions(solution, + explicit_solution, + solution_dot); + locally_relevant_solution = solution; + locally_relevant_explicit_solution = explicit_solution; + locally_relevant_explicit_solution_dot = solution_dot; + + + + // read zlib compressed resume.z + try + { +#ifdef DEAL_II_WITH_ZLIB + std::ifstream ifs ((snapshot_folder + "restart.resume.z.old").c_str()); + AssertThrow(ifs.is_open(), + ExcMessage("Cannot open snapshot resume file.")); + + uint32_t compression_header[4]; + ifs.read((char *)compression_header, 4 * sizeof(compression_header[0])); + Assert(compression_header[0]==1, ExcInternalError()); + + std::vector compressed(compression_header[3]); + std::vector uncompressed(compression_header[1]); + ifs.read(&compressed[0],compression_header[3]); + uLongf uncompressed_size = compression_header[1]; + + const int err = uncompress((Bytef *)&uncompressed[0], &uncompressed_size, + (Bytef *)&compressed[0], compression_header[3]); + AssertThrow (err == Z_OK, + ExcMessage (std::string("Uncompressing the data buffer resulted in an error with code <") + + + Utilities::int_to_string(err))); + + { + std::istringstream ss; + ss.str(std::string (&uncompressed[0], uncompressed_size)); + boost::archive::binary_iarchive ia (ss); + ia >> (*this); + } +#else + AssertThrow (false, + ExcMessage ("You need to have deal.II configured with the 'libz' " + "option to support checkpoint/restart, but deal.II " + "did not detect its presence when you called 'cmake'.")); +#endif + } + catch (std::exception &e) + { + AssertThrow (false, + ExcMessage (std::string("Cannot seem to deserialize the data previously stored!\n") + + + "Some part of the machinery generated an exception that says <" + + + e.what() + + + ">")); + } +} + + +template +void +piDoMUS::load_solutions(LATrilinos::VectorType &y, + LATrilinos::VectorType &y_expl, + LATrilinos::VectorType &y_dot) +{ + LATrilinos::VectorType distributed_system (y); + LATrilinos::VectorType expl_distributed_system (y_expl); + LATrilinos::VectorType distributed_system_dot (y_dot); + + std::vector x_system (3); + x_system[0] = & (distributed_system); + x_system[1] = & (expl_distributed_system); + x_system[2] = & (distributed_system_dot); + + parallel::distributed::SolutionTransfer > + system_trans (*dof_handler); + + system_trans.deserialize (x_system); + + y = distributed_system; + y_expl = expl_distributed_system; + y_dot = distributed_system_dot; + +} + + +template +void +piDoMUS::load_solutions(LADealII::VectorType &y, + LADealII::VectorType &y_expl, + LADealII::VectorType &y_dot) +{ + LADealII::VectorType distributed_system (y); + LADealII::VectorType expl_distributed_system (y_expl); + LADealII::VectorType distributed_system_dot (y_dot); + + std::vector x_system (3); + x_system[0] = distributed_system; + x_system[1] = expl_distributed_system; + x_system[2] = distributed_system_dot; + + SolutionTransfer > + system_trans (*dof_handler); + +//system_trans.deserialize (x_system); + + y = distributed_system; + y_expl = expl_distributed_system; + y_dot = distributed_system_dot; + +} + + template class piDoMUS<2, 2, LATrilinos>; template class piDoMUS<2, 3, LATrilinos>; template class piDoMUS<3, 3, LATrilinos>; diff --git a/tests/parameters/poisson_problem_01.prm b/tests/parameters/poisson_problem_01.prm index 61610798..6766cd03 100644 --- a/tests/parameters/poisson_problem_01.prm +++ b/tests/parameters/poisson_problem_01.prm @@ -1,8 +1,8 @@ # Parameter file generated with # D2K_GIT_BRANCH= master -# D2K_GIT_SHORTREV= bedec02 -# DEAL_II_GIT_BRANCH= pdsolution-transfer -# DEAL_II_GIT_SHORTREV= fe83377 +# D2K_GIT_SHORTREV= fc320d8 +# DEAL_II_GIT_BRANCH= master +# DEAL_II_GIT_SHORTREV= 372f54d subsection Dirichlet boundary conditions set IDs and component masks = 0=ALL set IDs and expressions = @@ -10,19 +10,23 @@ subsection Dirichlet boundary conditions set Used constants = end subsection Domain - set Colorize = false - set Grid to generate = rectangle - set Input grid file name = - set Mesh smoothing alogrithm = none - set Optional Point 1 = 1,2,7 - set Optional Point 2 = 3,4,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 = + 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 = 1,2,7 + set Optional Point 2 = 3,4,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 @@ -72,16 +76,19 @@ subsection IDA Solver Parameters end subsection IMEX Parameters set Absolute error tolerance = 1e-6 - set Final time = 0. + set Final time = 1. 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 Method used = fixed_alpha set Newton relaxation parameter = 1.000000 + set Number of elements in backtracking sequence = 5 + set Print useful informations = false set Relative error tolerance = 0.000000 - set Step size = 1e10 + set Step size = 1e-3 set Update continuously Jacobian = true - set Norm used for non linear iterations = l2 + set Use the KINSOL solver = true end subsection Initial solution set Function constants = @@ -93,6 +100,15 @@ subsection Initial solution_dot set Function expression = 0 set Variable names = x,y,z,t end +subsection KINSOL for IMEX + set Level of verbosity of the KINSOL solver = 0 + set Maximum number of iteration before Jacobian update = 10 + set Maximum number of iterations = 200 + set Step tolerance = 1e-11 + set Strategy = newton + set Tolerance for residuals = 1e-9 + set Use internal KINSOL direct solver = false +end subsection Neumann boundary conditions set IDs and component masks = set IDs and expressions = @@ -126,15 +142,27 @@ subsection Time derivative of Dirichlet boundary conditions 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 = 1 - set Maximum number of time steps = 10000 + set Enable finer preconditioner = false + set Initial global refinement = 3 set Jacobian solver tolerance = 1e-8 + set Max iterations = 50 + set Max iterations finer prec. = 0 + 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 Time stepper = euler + set Refine mesh during transient = false + set Resume computation from a snapshot = true set Show timer = false + set Snapshot folder = snapshots + set Threshold for solver's restart = 1e-2 + set Time stepper = euler set Use direct solver if available = true end