IDAInterface #4852
IDAInterface #4852
Conversation
|
@jdannberg here it is. :) |
luca-heltai
force-pushed
the
ida-interface
branch
from
05e48d7
to
2c5aae2
Compare
|
Great, thanks! |
| * scalar \f$\alpha\f$ changes whenever the step size or method order | ||
| * changes. | ||
| */ | ||
| template<typename VEC=Vector<double> > |
There was a problem hiding this comment.
Elsewhere we call template VectorType
include/deal.II/sundials/copy.h
Outdated
| #ifdef DEAL_II_WITH_MPI | ||
|
|
||
| #ifdef DEAL_II_WITH_TRILINOS | ||
| void copy(TrilinosWrappers::MPI::Vector &dst, const N_Vector &src); |
There was a problem hiding this comment.
This and below are not documented
| */ | ||
| IDAInterface(const std::string name="", | ||
| const MPI_Comm mpi_comm = MPI_COMM_WORLD); | ||
| /** House cleaning. */ |
There was a problem hiding this comment.
Space before and indentation the same way as constructor
| /** House cleaning. */ | ||
| ~IDAInterface(); | ||
|
|
||
| /** Declare parameters for this class to function properly. */ |
There was a problem hiding this comment.
I don't understand what u mean by "to function properly"? Can u please reprhrase?
There was a problem hiding this comment.
looking at the test, i would list here explicitly all the parameters which will be declared. Otherwise users are forced to dig into the source code to see what is actually happening inside to be able to use this function.
I would take the output of ParameterHandler::print_parameters().
| const unsigned int step_number)> output_step; | ||
|
|
||
| /** | ||
| * Evaluate wether the solver should be restarted (for example because the |
| */ | ||
| void set_functions_to_trigger_an_assert(); | ||
|
|
||
| /** Initial time for the DAE.*/ |
There was a problem hiding this comment.
Could u please indent it the way we do elsewhere, i.e.
/**
* Initial time for the DAE.
*/
Same below.
| N_Vector diff_id; | ||
|
|
||
| #ifdef DEAL_II_WITH_MPI | ||
| MPI_Comm communicator; |
| /** Output stream */ | ||
| ConditionalOStream pcout; | ||
|
|
||
| unsigned int system_size; |
|
|
||
| unsigned int system_size; | ||
|
|
||
| unsigned int local_system_size; |
|
@jdannberg -- want to practice reviewing patches with this one? ;-) |
| /** House cleaning. */ | ||
| ~IDAInterface(); | ||
|
|
||
| /** Declare parameters for this class to function properly. */ |
There was a problem hiding this comment.
looking at the test, i would list here explicitly all the parameters which will be declared. Otherwise users are forced to dig into the source code to see what is actually happening inside to be able to use this function.
I would take the output of ParameterHandler::print_parameters().
source/sundials/ida_interface.cc
Outdated
| { | ||
| prm.add_parameter("Initial step size",initial_step_size); | ||
|
|
||
| prm.add_parameter("Min step size", min_step_size); |
There was a problem hiding this comment.
since you say Maximum elsewhere, please use Minimum here.
source/sundials/ida_interface.cc
Outdated
|
|
||
| prm.add_parameter("Final time", final_time); | ||
|
|
||
| prm.add_parameter("Time units between each output", output_period); |
There was a problem hiding this comment.
not being a native speaker, i perceive units as units of time (second, hours, etc), thus the description is confusing IMO. Maybe change to Time interval between each output?
source/sundials/ida_interface.cc
Outdated
| "in the local error test."); | ||
|
|
||
| prm.add_parameter("Initial condition type", ic_type, | ||
| "This is one of the following thress options for the " |
source/sundials/ida_interface.cc
Outdated
| Patterns::Selection("none|use_y_diff|use_y_dot")); | ||
|
|
||
| prm.add_parameter("Initial condition type after restart", reset_type, | ||
| "This is one of the following thress options for the " |
|
@luca-heltai please also add a quick test for the new interface, see #3230 |
|
I'm still off the grid with only my phone (and my family :-)). As soon as I'm back to Trieste I'll fix all comments, and add some more stuff. |
|
OK, no rush. Enjoy the time! |
|
I think I have addressed all of @davydden requests. If this is looking good, I'll add to sundials directory also the interfaces to the general nonlinear solver ( |
include/deal.II/sundials/copy.h
Outdated
| @@ -0,0 +1,74 @@ | |||
| //----------------------------------------------------------- | |||
| // | |||
| // Copyright (C) 2015 by the deal.II authors | |||
There was a problem hiding this comment.
You may want to adjust the year.
Do we want to stick with the directory name sundials/? Or should we go with something more generic, say odesolvers/?
There was a problem hiding this comment.
sundials offers also nonlinear solvers and dae solvers (IDA is a dae solver, for example, not an ode solver...). I'd stick to sundials.
There was a problem hiding this comment.
i agree with @luca-heltai . Eventually we may have more interfaces to other libs which contain different types of solvers so i would group headers based on the library name.
| #ifdef DEAL_II_WITH_MPI | ||
|
|
||
| #ifdef DEAL_II_WITH_TRILINOS | ||
| void copy(TrilinosWrappers::MPI::Vector &dst, const N_Vector &src); |
There was a problem hiding this comment.
You choose the vmult convention of saying copy(dst,src) (which mirrors dst = src). But that's not what std::copy does, for example, nor what most of the rest of the deal.II functions do. Thoughts?
There was a problem hiding this comment.
I kept this in an internal namespace, because ideally we should get rid of all these in favour of ReadWriteVector and VectorSpaceVector objects. All copy functions should disappear.
| * \end{cases} | ||
| * \f] | ||
| * | ||
| * where \f$y,\dot y\f$ are vectors in \f$\R^n\f$, \f$t\f$ is often the time (but can |
There was a problem hiding this comment.
You don't need to use \f (in fact, it is wrong) since we automatically expand $ to \f$.
There was a problem hiding this comment.
I copyed from deal2lkit... there we don't have this auto-magic expansion... :) I'll fix this.
| * also be a parametric quantity), and | ||
| * \f$F:\R\times\R^n\times\R^n\rightarrow\R^n\f$. Such problem is solved | ||
| * using Newton iteration augmented with a line search global | ||
| * strategy. The integration method used in ida is the variable-order, |
| * the BDF of order \f$q\f$ given by the multistep formula | ||
| * | ||
| * \f[ | ||
| * \sum\limits_{i=0}^q \alpha_{n,i}\,y_{n-i}=h_n\,\dot y_n\, , |
There was a problem hiding this comment.
\limit? (or maybe I just don't know the command?)
source/sundials/ida_interface.cc
Outdated
| use_local_tolerances(use_local_tolerances), | ||
| ida_mem(nullptr), | ||
| communicator(Utilities::MPI::duplicate_communicator(mpi_comm)), | ||
| pcout(std::cout, Utilities::MPI::this_mpi_process(mpi_comm)==0) |
There was a problem hiding this comment.
Could we do this differently, for example via a signal?
source/sundials/ida_interface.cc
Outdated
| next_time += output_period; | ||
| if (verbose) | ||
| { | ||
| pcout << " "//"\r" |
source/sundials/ida_interface.cc
Outdated
| << std::setw(5) << next_time | ||
| << std::endl; | ||
| } | ||
| status = IDASolve(ida_mem, next_time, &t, yy, yp, IDA_NORMAL); |
There was a problem hiding this comment.
can you add a check of the return value here?
source/sundials/ida_interface.cc
Outdated
| status = IDASolve(ida_mem, next_time, &t, yy, yp, IDA_NORMAL); | ||
|
|
||
| status = IDAGetLastStep(ida_mem, &h); | ||
| AssertThrow(status == 0, ExcMessage("Error in IDA Solver")); |
There was a problem hiding this comment.
I think it would be nice to output the error code. I think there is already something like AssertMPI or Assert(PETSc) that you could steal for here.
source/sundials/ida_interface.cc
Outdated
| // double frac = 0; | ||
| int k = 0; | ||
| IDAGetLastOrder(ida_mem, &k); | ||
| // frac = std::pow((double)k,2.); |
luca-heltai
force-pushed
the
ida-interface
branch
2 times, most recently
from
16ab532
to
d9cf16c
Compare
|
I think this is ready for review. As soon as I have some more time, I'll add more of the functionality provided in IDA to this class (with additional tests). |
include/deal.II/sundials/ida.h
Outdated
| { | ||
|
|
||
| /** | ||
| * Interface to SUNDIALS IDA (Implicit Differential-Algebraic) solver. |
There was a problem hiding this comment.
maybe introduce abbreviation here:Implicit Differential-Algebraic (IDA) and then you can use it below?
include/deal.II/sundials/ida.h
Outdated
| const bool &use_local_tolerances = false); | ||
|
|
||
| /** | ||
| * House cleaning. |
There was a problem hiding this comment.
for consistency with other classes, i would write Destructor.
include/deal.II/sundials/ida.h
Outdated
| virtual void add_parameters(ParameterHandler &prm); | ||
|
|
||
| /** | ||
| * Evolve. This function returns the final number of steps. |
There was a problem hiding this comment.
maybe Evolve -> Integrate differential-dlgebraic equations?
include/deal.II/sundials/ida.h
Outdated
| * called when the simulation start and when the user returns true to a | ||
| * call to solver_should_restart(). | ||
| * | ||
| * By default solver_should_restart returns false. If the user needs to |
include/deal.II/sundials/ida.h
Outdated
| /** | ||
| * Minimum step size. | ||
| */ | ||
| double min_step_size; |
There was a problem hiding this comment.
for consistency with the two above minimum_step_size?
include/deal.II/sundials/ida.h
Outdated
| /** | ||
| * Absolute error tolerance for adaptive time stepping. | ||
| */ | ||
| double abs_tol; |
include/deal.II/sundials/ida.h
Outdated
| /** | ||
| * Relative error tolerance for adaptive time stepping. | ||
| */ | ||
| double rel_tol; |
include/deal.II/sundials/ida.h
Outdated
| /** | ||
| * Maximum order of BDF. | ||
| */ | ||
| unsigned int max_order; |
include/deal.II/sundials/ida.h
Outdated
| const unsigned int &max_order = 5, | ||
| const double &output_period = .1, | ||
| const bool &ignore_algebraic_terms_for_errors = true, | ||
| const std::string &ic_type = "none", |
There was a problem hiding this comment.
i would switch ic_type and reset_typeto enumerators and only internally convert to std::string to provide the correct input to Sundials. Same for class private variables.
include/deal.II/sundials/ida.h
Outdated
| * Output stream. If run in parallel, only processor zero will | ||
| * output verbose information about time steps. | ||
| */ | ||
| ConditionalOStream pcout; |
There was a problem hiding this comment.
@dealii/dealii if I understand correctly we are moving away form writing anything from dealii's classes and instead provide functions for optional output
include/deal.II/sundials/ida.h
Outdated
| * Those options are set as default value in the ParameterHandler object | ||
| * `prm`. | ||
| */ | ||
| virtual void add_parameters(ParameterHandler &prm); |
There was a problem hiding this comment.
an alternative is to move those into an auxiliary class AdditionalData which can declare and read things from ParameterHandler and has a constructor with default values. Then you can easily extend the number of parameters later without breaking constructor interface of the main class.
What I find more important is that you can then make all private fixed parameter members const.
include/deal.II/sundials/ida.h
Outdated
| * | ||
| * That is $F(y', y, t) = y' + A y = 0 $ | ||
| * | ||
| * where |
include/deal.II/sundials/ida.h
Outdated
| * \end{matrix} | ||
| * \f] | ||
| * | ||
| * and $y(0)=(0, k)$, $y'(0) = (k, 0)$ |
There was a problem hiding this comment.
same here.
period at the end of the sentence
include/deal.II/sundials/ida.h
Outdated
| * and $y(0)=(0, k)$, $y'(0) = (k, 0)$ | ||
| * | ||
| * The exact solution is $y_0(t) = \sin(k t)$, $y_1(t) = y_0'(t) = k \cos(k t)$, | ||
| * $y_1'(t) = -k^2 \sin(k t)$ |
There was a problem hiding this comment.
period. same in the following line
| std::function<int(const double t, | ||
| const VectorType &y, | ||
| const VectorType &y_dot, | ||
| const double alpha)> setup_jacobian; |
There was a problem hiding this comment.
since the function only takes const arguments, how does it actually setup a jacobian? Can you say when this function is called, and where the jacobian it sets up is stored?
There was a problem hiding this comment.
it is the user's responsability to decide where to store the Jacobian. This class is agnostic on what the user wants to do, and does not enforce it. I'll document this a little better. The user should make sure that they have a way of solving for J = alpha dF/dy_dot + dF/dy after a call to this function is made. If they want to store a matrix, or setup a matrix free computation, it is their responsability, and this class should not care...
There was a problem hiding this comment.
This is now documented.
| * - <0: Unrecoverable error the computation will be aborted and an assertion | ||
| * will be thrown. | ||
| */ | ||
| std::function<int(const VectorType &rhs, VectorType &dst)> solve_jacobian_system; |
There was a problem hiding this comment.
can you say what the function arguments mean, and that the function is supposed to compute dst = J^{-1} rhs?
There was a problem hiding this comment.
done. I don't know why it does not show up here.
| std::function<void (const double t, | ||
| const VectorType &sol, | ||
| const VectorType &sol_dot, | ||
| const unsigned int step_number)> output_step; |
There was a problem hiding this comment.
I think you wanted to update the description (and maybe name) to state when the function will be called and what it is supposed to do.
luca-heltai
force-pushed
the
ida-interface
branch
from
ef77868
to
0db9147
Compare
|
I've rebased and squashed to a smaller number of commits. Not as good as @drwells git-history (from his commits it looks like he always get things right the first time! :) ...) but this should be small enough. I think we are now ready. |
|
OK to merge once the tester is happy. /run-tests |
|
In reference to #5027. |
#4704