Merge pull request #207 from asartori86/no-sundials

piDoMUS is no longer derived from SundialsInterface

include/base_interface.h

@@ -399,7 +399,11 @@ class BaseInterface : public ParameterAcceptor, public SimulatorAccess<dim,space
   */
   void set_stepper(const std::string &s) const;
 
-
+  /**
+    * Return norm of vector @p v. This function is called by the IMEX stepper.
+    * By default it returns the l2 norm.
+    */
+  virtual double vector_norm(const typename LAC::VectorType &v) const;
 
 protected:
 

include/lac/lac_type.h

@@ -17,7 +17,7 @@
 
 
 #include <deal.II/base/config.h>
-
+#include <deal2lkit/config.h>
 
 #include <deal.II/lac/vector.h>
 #include <deal.II/lac/block_vector.h>
@@ -28,7 +28,7 @@
 #include <deal.II/lac/precondition.h>
 
 using namespace dealii;
-using namespace deal2lkit;
+//using namespace deal2lkit;
 
 class LADealII
 {

include/pidomus.h

@@ -34,14 +34,15 @@
 #include "base_interface.h"
 #include "simulator_access.h"
 #include "pidomus_signals.h"
+#include "pidomus_lambdas.h"
+
 
 #include <deal2lkit/parsed_grid_generator.h>
 #include <deal2lkit/parsed_finite_element.h>
 #include <deal2lkit/parsed_grid_refinement.h>
 #include <deal2lkit/error_handler.h>
 #include <deal2lkit/parsed_function.h>
 #include <deal2lkit/parameter_acceptor.h>
-#include <deal2lkit/sundials_interface.h>
 #include <deal2lkit/ida_interface.h>
 #include <deal2lkit/imex_stepper.h>
 #include <deal2lkit/parsed_zero_average_constraints.h>
@@ -59,7 +60,7 @@ using namespace deal2lkit;
 using namespace pidomus;
 
 template <int dim, int spacedim = dim, typename LAC = LATrilinos>
-class piDoMUS : public ParameterAcceptor, public SundialsInterface<typename LAC::VectorType>
+class piDoMUS : public ParameterAcceptor
 {
 
 
@@ -69,10 +70,14 @@ class piDoMUS : public ParameterAcceptor, public SundialsInterface<typename LAC:
 
 public:
 
-
+#ifdef DEAL_II_WITH_MPI
   piDoMUS (const std::string &name,
            const BaseInterface<dim, spacedim, LAC> &energy,
-           const MPI_Comm &comm = MPI_COMM_WORLD);
+           const MPI_Comm comm = MPI_COMM_WORLD);
+#else
+  piDoMUS (const std::string &name,
+           const BaseInterface<dim, spacedim, LAC> &energy);
+#endif
 
   virtual void declare_parameters(ParameterHandler &prm);
   virtual void parse_parameters_call_back();
@@ -96,17 +101,14 @@ class piDoMUS : public ParameterAcceptor, public SundialsInterface<typename LAC:
   virtual void output_step(const double t,
                            const typename LAC::VectorType &solution,
                            const typename LAC::VectorType &solution_dot,
-                           const unsigned int step_number,
-                           const double h);
+                           const unsigned int step_number);
 
   /** 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,
                                      typename LAC::VectorType &solution,
                                      typename LAC::VectorType &solution_dot);
 
@@ -121,17 +123,11 @@ class piDoMUS : public ParameterAcceptor, public SundialsInterface<typename LAC:
   virtual int setup_jacobian(const double t,
                              const typename LAC::VectorType &src_yy,
                              const typename LAC::VectorType &src_yp,
-                             const typename LAC::VectorType &residual,
                              const double alpha);
 
 
   /** Inverse of the Jacobian vector product. */
-  virtual int solve_jacobian_system(const double t,
-                                    const typename LAC::VectorType &y,
-                                    const typename LAC::VectorType &y_dot,
-                                    const typename LAC::VectorType &residual,
-                                    const double alpha,
-                                    const typename LAC::VectorType &src,
+  virtual int solve_jacobian_system(const typename LAC::VectorType &src,
                                     typename LAC::VectorType &dst) const;
 
   /**
@@ -307,13 +303,14 @@ class piDoMUS : public ParameterAcceptor, public SundialsInterface<typename LAC:
 
   void set_constrained_dofs_to_zero(typename LAC::VectorType &v) const;
 
-  const MPI_Comm &comm;
+#ifdef DEAL_II_WITH_MPI
+  const MPI_Comm comm;
+#endif
   const BaseInterface<dim, spacedim, LAC>    &interface;
 
   unsigned int n_cycles;
   unsigned int current_cycle;
   unsigned int initial_global_refinement;
-  unsigned int max_time_iterations;
 
   ConditionalOStream        pcout;
 
@@ -436,7 +433,7 @@ class piDoMUS : public ParameterAcceptor, public SundialsInterface<typename LAC:
 
 
   IDAInterface<typename LAC::VectorType>  ida;
-  IMEXStepper<typename LAC::VectorType> euler;
+  IMEXStepper<typename LAC::VectorType> imex;
 
 
   std::vector<types::global_dof_index> dofs_per_block;
@@ -522,6 +519,12 @@ class piDoMUS : public ParameterAcceptor, public SundialsInterface<typename LAC:
                   const typename LAC::VectorType &solution_dot);
 
 
+  /**
+   * Return the norm of vector @p v. Internally it calls the
+   * virtual function BaseInterface::vector_norm().
+   */
+  double vector_norm(const typename LAC::VectorType &v) const;
+
   /**
    * Struct containing the signals
    */
@@ -533,6 +536,12 @@ class piDoMUS : public ParameterAcceptor, public SundialsInterface<typename LAC:
    */
   friend class SimulatorAccess<dim,spacedim,LAC>;
 
+public:
+
+  friend class Lambdas<dim,spacedim,LAC>;
+
+  Lambdas<dim,spacedim,LAC> lambdas;
+
 };
 
 

include/pidomus_lambdas.h

@@ -0,0 +1,148 @@
+#ifndef __pidomus_lambdas_h_
+#define __pidomus_lambdas_h_
+
+#include "lac/lac_type.h"
+#include <deal2lkit/utilities.h>
+
+using namespace deal2lkit;
+// forward declaration
+template <int dim, int spacedim, typename LAC> class piDoMUS;
+
+/**
+ * Lambdas class. A piDoMUS object has a Lambdas object and the
+ * std::functions of the provided stepper (e.g., ida, imex)
+ * are set to be equal to those here implemented.
+ *
+ * By default, the std::functions of this class call the
+ * namesake functions of the piDoMUS object, which is specified
+ * either when the Lambdas object is constructed or by calling the
+ * Lambdas::initialize_simulator() function.
+ *
+ * The aim of this class is to increase the flexibility of piDoMUS.
+ * piDoMUS offers flexibility by itself thanks to the signals
+ * to whom the user can connect in order to perform problem-specific
+ * tasks. Whether the behavior of a function should be completely different,
+ * the user can ovveride the functions declared here (without the need
+ * of modifying piDoMUS' source code).
+ */
+
+template <int dim, int spacedim, typename LAC>
+class Lambdas
+{
+public:
+
+  /**
+   * Default constructor. Initialize the Lambdas object without
+   * a reference to a particular piDoMUS object. You will later have
+   * to call initialize() to provide this reference to the piDoMUS
+   * object.
+   */
+  Lambdas ();
+
+  /**
+   * Create a Lambdas object that is already initialized for
+   * a particular piDoMUS.
+   */
+  Lambdas (piDoMUS<dim,spacedim,LAC> &simulator_object);
+
+  /**
+   * Destructor. Does nothing.
+   */
+  ~Lambdas ();
+
+  /**
+   * Initialize this class for a given simulator.
+   *
+   * @param simulator_object A reference to the main simulator object.
+   */
+  void initialize_simulator (piDoMUS<dim, spacedim, LAC> &simulator_object);
+
+  /**
+   * Set the default behavior of the functions of this class,
+   * which is call the namesake functions of the simulator_bject.
+   */
+  void set_functions_to_default();
+
+  /**
+   * Return a shared_ptr<VECTOR_TYPE>. A shared_ptr is needed in order to
+   * keep the pointed vector alive, without the need to use a static variable.
+   */
+  std::function<shared_ptr<typename LAC::VectorType>()> create_new_vector;
+
+  /**
+   * Compute residual.
+   */
+  std::function<int(const double t,
+                    const typename LAC::VectorType &y,
+                    const typename LAC::VectorType &y_dot,
+                    typename LAC::VectorType &res)> residual;
+
+  /**
+   * Compute Jacobian.
+   */
+  std::function<int(const double t,
+                    const typename LAC::VectorType &y,
+                    const typename LAC::VectorType &y_dot,
+                    const double alpha)> setup_jacobian;
+
+  /** Solve linear system. */
+  std::function<int(const typename LAC::VectorType &rhs, typename LAC::VectorType &dst)> solve_jacobian_system;
+
+  /**
+   * Store solutions to file.
+   */
+  std::function<void (const double t,
+                      const typename LAC::VectorType &sol,
+                      const typename LAC::VectorType &sol_dot,
+                      const unsigned int step_number)> output_step;
+
+  /**
+   * Evaluate wether the mesh should be refined or not. If so,
+   * it refines and interpolate the solutions from the old to the
+   * new mesh.
+   */
+  std::function<bool (const double t,
+                      typename LAC::VectorType &sol,
+                      typename LAC::VectorType &sol_dot)> solver_should_restart;
+
+  /**
+   * Return a vector whose component are 1 if the corresponding
+   * dof is differential, 0 if algebraic. This function is needed
+   * by the IDAInterface stepper.
+   */
+  std::function<typename LAC::VectorType&()> differential_components;
+
+  /**
+   * Return a vector whose components are the weights used by
+   * IDA to compute the vector norm. By default this function is not
+   * implemented.
+   */
+  std::function<typename LAC::VectorType&()> get_local_tolerances;
+
+  /**
+   * Return a vector which is a lumped mass matrix. This function
+   * is used by Kinsol (through imex) for setting the weights used
+   * for computing the norm a vector.
+   */
+  std::function<typename LAC::VectorType&()> get_lumped_mass_matrix;
+
+  /**
+   * Compute the matrix-vector product Jacobian times @p src,
+   * and the result is put in @p dst. It is required by Kinsol.
+   */
+  std::function<int(const typename LAC::VectorType &src,
+                    typename LAC::VectorType &dst)> jacobian_vmult;
+
+  /**
+   * Return the norm of @p vector, which is used by IMEXStepper.
+   */
+  std::function<double(const typename LAC::VectorType &vector)> vector_norm;
+
+  /**
+   * A pointer to the simulator object to which we want to get
+   * access.
+   */
+  piDoMUS<dim,spacedim,LAC> *simulator;
+};
+
+#endif

include/simulator_access.h

@@ -31,7 +31,7 @@ class SimulatorAccess
   SimulatorAccess ();
 
   /**
-   * Create a piDoMUSAccess object that is already initialized for
+   * Create a SimulatorAccess object that is already initialized for
    * a particular piDoMUS.
    */
   SimulatorAccess (const piDoMUS<dim,spacedim,LAC> &simulator_object);
@@ -73,11 +73,13 @@ class SimulatorAccess
   Signals<dim,spacedim,LAC> &
   get_signals() const;
 
+#ifdef DEAL_II_WITH_MPI
   /**
    * Return the MPI communicator for this simulation.
    */
   MPI_Comm
   get_mpi_communicator () const;
+#endif
 
   /**
    * Return a reference to the stream object that only outputs something

source/base_interface.cc

@@ -552,6 +552,15 @@ void
 BaseInterface<dim,spacedim,LAC>::connect_to_signals() const
 {}
 
+template<int dim, int spacedim, typename LAC>
+double
+BaseInterface<dim,spacedim,LAC>::
+vector_norm (const typename LAC::VectorType &v) const
+{
+  return v.l2_norm();
+}
+
+
 #define INSTANTIATE(dim,spacedim,LAC) \
   template class BaseInterface<dim,spacedim,LAC>;