Merge pull request #125 from Pressio/develop

v0.14.0
Pressio · Jun 14, 2023 · 682625a · 682625a
2 parents 312033a + ba9340e
commit 682625a
Show file tree

Hide file tree

Showing 306 changed files with 11,212 additions and 14,278 deletions.
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -3,17 +3,26 @@ cmake_minimum_required(VERSION 3.18.0)
 project(pressiodemoapps CXX)
 
 #=====================================================================
-# we need c++14
-set(CMAKE_CXX_STANDARD 14)
+# versioning
+#=====================================================================
+file(READ "version.txt" PDA_VERSION)
+message("pressiodemoapps version = ${PDA_VERSION}")
+
+#=====================================================================
+# we need c++17
+set(CMAKE_CXX_STANDARD 17)
 set(CMAKE_CXX_STANDARD_REQUIRED ON)
 set(CMAKE_CXX_EXTENSIONS OFF)
 # test if compiler supports standard
 include(CheckCXXCompilerFlag)
-check_cxx_compiler_flag("-std=c++14" COMPILER_SUPPORT_CPP14)
-if(NOT COMPILER_SUPPORT_CPP14)
-  message(FATAL_ERROR "Compiler does not support -std=c++14. This is required.")
+check_cxx_compiler_flag("-std=c++17" COMPILER_SUPPORT_CPP17)
+if(NOT COMPILER_SUPPORT_CPP17)
+  message(FATAL_ERROR "Compiler does not support -std=c++17. This is required.")
 endif()
 
+# force C++17 for pressio too
+set(PRESSIO_ENABLE_CXX17 ON)
+
 #=====================================================================
 # guard against in-source builds
 if(${CMAKE_SOURCE_DIR} STREQUAL ${CMAKE_BINARY_DIR})

diff --git a/README.md b/README.md
@@ -4,13 +4,14 @@
 **pressio-demoapps** is a library of 1D, 2D and 3D demo problems
 of varying complexity, ranging from a simple 1D linear advection,
 to 2D reaction-diffusion, and 3D Euler equations, and more.
+
 Key features include:
 
 - support for both C++ and Python
 - cell-centered finite volume discretization with various numerical schemes and *exact Jacobians*
 - focus on providing self-contained and well-defined problems
 - built-in support for a sample mesh: this mean that one can evaluate the residual and Jacobian
-at a disjoint subset of the mesh cells (this is useful for intrusive ROMs)
+at a disjoint subset of the mesh cells (this is useful for intrusive ROMs, but in general for other purposes too)
 
 Click below to check the documentation for more details:
 
@@ -20,10 +21,9 @@ Click below to check the documentation for more details:
 
 ## Development status
 
-**pressio-demoapps** is planned to be a long-term maintained project.
-Therefore, more problems and features will be implemented.
-If you are interested in collaborating or would like to see a specific
-problem added, please reach out.
+**pressio-demoapps** is planned to be maintained in the long-term.
+Hopefully, more problems and features will be implemented.
+If you are interested in collaborating or would like to see a specific problem added, please reach out.
 
 ## Questions?
 Find us on Slack: https://pressioteam.slack.com or open an issue on [github](https://github.com/Pressio/pressio-demoapps).

diff --git a/docs/source/apicpp.rst b/docs/source/apicpp.rst
@@ -38,6 +38,10 @@ A *pressio-demoapps* C++ problem class meets the following API
       (see `Supported Backends And Types`_ for more details about backends and corresponding types).
       For example, if you use ``create_problem_eigen``, the ``right_hand_side_type`` will be an Eigen vector.
 
+   .. cpp:type:: rhs_type
+
+      Shortcut for ``right_hand_side_type``
+
    .. cpp:type:: jacobian_type
 
       Data structure type to store the Jacobian: this depends on which specific implementation
@@ -53,6 +57,10 @@ A *pressio-demoapps* C++ problem class meets the following API
 
       Constructs and returns an instance of the right hand side.
 
+   .. cpp:function:: right_hand_side_type createRhs()
+
+      Shortcut for ``createRightHandSide()``.
+
    .. cpp:function:: jacobian_type createJacobian()
 
       Constructs and returns an instance of the jacobian.
@@ -69,6 +77,13 @@ A *pressio-demoapps* C++ problem class meets the following API
       evaluates the RHS of the system overwriting :math:`v` and  if ``computeJac == true``,
       its Jacobian and stores it into :math:`J`.
 
+   .. cpp:function:: void rhsAndJacobian(const state_type & y, scalar_type time,\
+       rhs_type & v, std::optional<jacobian_type *> J)
+
+      Given a state :math:`y` and time :math:`time`,
+      evaluates the RHS of the system overwriting :math:`v` and if ``J == true``,
+      its Jacobian and stores it into ``*(J.value())``.
+
    .. cpp:function:: auto totalDofSampleMesh()
 
       Returns the total number of degrees of freedom on the **sample** mesh.

diff --git a/docs/source/installation.rst b/docs/source/installation.rst
@@ -1,15 +1,14 @@
 .. highlight:: sh
 
-Installation
-============
-
+Build/install
+=============
 
 C++ library
 -----------
 
-The C++ library is header-only so it does not need to be compiled and installed.
+**The C++ library is header-only so it does not need to be compiled and installed.**
 
-To use it, you need a C++14 compiler and you have to:
+To use it, you need a C++17 compiler and you have to:
 
 1. include the ``pressiodemoapps/include`` subdirectory in your compilation line
 
@@ -21,7 +20,7 @@ Building the test suite
 If you want to build the C++ tests, you need CMake > 3.18.0 and then do::
 
   git clone --recursive [email protected]:Pressio/pressio-demoapps.git
-  export CXX=<path-to-your-CXX-compiler> #must support C++14
+  export CXX=<path-to-your-CXX-compiler> #must support C++17
   cd pressio-demoapps && mkdir build && cd build
   cmake -DPRESSIODEMOAPPS_ENABLE_TESTS=On ..
   make -j4
@@ -33,7 +32,7 @@ Python library
 
 To build/install the bindings, you need:
 
-- CMake > 3.18.0 and a C++ compiler with C++14 support
+- CMake > 3.18.0 and a C++ compiler with C++17 support
 
 - You also need these packages::
 
@@ -44,7 +43,7 @@ Then, you can do:
 .. code-block:: shell
 
    git clone --recursive [email protected]:Pressio/pressio-demoapps.git
-   export CXX=<path-to-your-CXX-compiler> #must support C++14
+   export CXX=<path-to-your-CXX-compiler> #must support C++17
    cd pressio-demoapps
    python3 cmake_build.py
    pip3 install .

diff --git a/include/pressiodemoapps/adapter_cpp.hpp b/include/pressiodemoapps/adapter_cpp.hpp
@@ -49,6 +49,8 @@
 #ifndef PRESSIODEMOAPPS_ADAPTER_MIXINS_CPP_HPP_
 #define PRESSIODEMOAPPS_ADAPTER_MIXINS_CPP_HPP_
 
+#include <optional>
+
 namespace pressiodemoapps{
 
 // this is the class that gets instantiated and passed
@@ -63,6 +65,7 @@ class PublicProblemEigenMixinCpp : public T
   using independent_variable_type = typename T::scalar_type;
   using time_type                 = typename T::scalar_type;
   using state_type		  = typename T::state_type;
+  using rhs_type                  = typename T::velocity_type;
   using right_hand_side_type      = typename T::velocity_type;
   using jacobian_type             = typename T::jacobian_type;
 
@@ -102,6 +105,11 @@ class PublicProblemEigenMixinCpp : public T
     return pda::clone(m_rhs);
   }
 
+  rhs_type createRhs() const{
+    namespace pda = ::pressiodemoapps;
+    return pda::clone(m_rhs);
+  }
+
   jacobian_type createJacobian() const{
     namespace pda = ::pressiodemoapps;
     return pda::clone(m_jacobian);
@@ -166,6 +174,14 @@ class PublicProblemEigenMixinCpp : public T
     }
   }
 
+  void rhs(const state_type & state,
+	   const independent_variable_type currentTime,
+	   right_hand_side_type & V) const
+  {
+    T::velocityAndOptionalJacobian(state, currentTime,
+				   V, nullptr);
+  }
+
   void rightHandSide(const state_type & state,
 		     const independent_variable_type currentTime,
 		     right_hand_side_type & V) const
@@ -183,6 +199,19 @@ class PublicProblemEigenMixinCpp : public T
 				   V, &jacobian);
   }
 
+  void rhsAndJacobian(const state_type & state,
+		      const independent_variable_type currentTime,
+		      rhs_type & V,
+		      std::optional<jacobian_type*> jacobian) const
+  {
+    if (jacobian){
+      T::velocityAndOptionalJacobian(state, currentTime, V, jacobian.value());
+    }
+    else{
+      T::velocityAndOptionalJacobian(state, currentTime, V, nullptr);
+    }
+  }
+
   void jacobian(const state_type & state,
 		const independent_variable_type currentTime,
 		jacobian_type & jacobian) const

diff --git a/meshing_scripts/create_full_mesh_for.py b/meshing_scripts/create_full_mesh_for.py
@@ -64,9 +64,6 @@
          "swe2dSlipWall_s3",\
          "swe2dSlipWall_s5",\
          "swe2dSlipWall_s7",\
-         "swe2dslipwall_s3",\
-         "swe2dslipwall_s5",\
-         "swe2dslipwall_s7",\
          #
          "burgers2d_periodic_s3", \
          "burgers2d_periodic_s5", \

diff --git a/tests_cpp/eigen_1d_diffusion_reaction_implicit/main.cc b/tests_cpp/eigen_1d_diffusion_reaction_implicit/main.cc
@@ -20,9 +20,8 @@ int main()
   using lin_solver_t = pressio::linearsolvers::Solver<
     pressio::linearsolvers::iterative::Bicgstab, jacob_t>;
   lin_solver_t linSolverObj;
-  auto NonLinSolver=
-   pressio::nonlinearsolvers::create_newton_raphson(stepperObj, linSolverObj);
-  NonLinSolver.setTolerance(1e-11);
+  auto NonLinSolver=pressio::create_newton_solver(stepperObj, linSolverObj);
+  NonLinSolver.setStopTolerance(1e-11);
 
   FomObserver<state_t> Obs("1d_diffreac_solution.bin", 10);
 

diff --git a/tests_cpp/eigen_1d_euler_sod_implicit/main.cc b/tests_cpp/eigen_1d_euler_sod_implicit/main.cc
@@ -34,9 +34,8 @@ int main()
     pressio::linearsolvers::iterative::Bicgstab, jacob_t>;
   lin_solver_t linSolverObj;
 
-  auto NonLinSolver=
-    pressio::nonlinearsolvers::create_newton_raphson(stepperObj, linSolverObj);
-  NonLinSolver.setTolerance(1e-6);
+  auto NonLinSolver=pressio::create_newton_solver(stepperObj, linSolverObj);
+  NonLinSolver.setStopTolerance(1e-6);
 
   FomObserver<state_t> Obs("sod1d_solution.bin", 1);
 

diff --git a/tests_cpp/eigen_1d_linear_advection_default_velocity_implicit/main.cc b/tests_cpp/eigen_1d_linear_advection_default_velocity_implicit/main.cc
@@ -36,9 +36,8 @@ int main()
     pressio::linearsolvers::iterative::Bicgstab, jacob_t>;
   lin_solver_t linSolverObj;
 
-  auto NonLinSolver=
-    pressio::nonlinearsolvers::create_newton_raphson(stepperObj, linSolverObj);
-  NonLinSolver.setTolerance(1e-6);
+  auto NonLinSolver= pressio::create_newton_solver(stepperObj, linSolverObj);
+  NonLinSolver.setStopTolerance(1e-6);
 
   FomObserver<state_t> Obs("linadv_1d_solution.bin", 1);
 

diff --git a/tests_cpp/eigen_2d_euler_double_mach_reflection_implicit/main.cc b/tests_cpp/eigen_2d_euler_double_mach_reflection_implicit/main.cc
@@ -34,8 +34,8 @@ int main()
   using lin_solver_t = pressio::linearsolvers::Solver<
     pressio::linearsolvers::iterative::Bicgstab, jacob_t>;
   lin_solver_t linSolverObj;
-  auto NonLinSolver = pressio::nonlinearsolvers::create_newton_raphson(stepperObj, linSolverObj);
-  NonLinSolver.setTolerance(1e-5);
+  auto NonLinSolver = pressio::create_newton_solver(stepperObj, linSolverObj);
+  NonLinSolver.setStopTolerance(1e-5);
 
   FomObserver<state_t> Obs("doubleMach2d_solution.bin", 10);
 

diff --git a/tests_cpp/eigen_2d_euler_normal_shock_implicit/main.cc b/tests_cpp/eigen_2d_euler_normal_shock_implicit/main.cc
@@ -34,8 +34,8 @@ int main()
   using lin_solver_t = pressio::linearsolvers::Solver<
     pressio::linearsolvers::iterative::Bicgstab, jacob_t>;
   lin_solver_t linSolverObj;
-  auto NonLinSolver= pressio::nonlinearsolvers::create_newton_raphson(stepperObj, linSolverObj);
-  NonLinSolver.setTolerance(1e-5);
+  auto NonLinSolver= pressio::create_newton_solver(stepperObj, linSolverObj);
+  NonLinSolver.setStopTolerance(1e-5);
 
   FomObserver<state_t> Obs("normalShock2d_solution.bin", 50);
 

diff --git a/tests_cpp/eigen_2d_euler_sedov_implicit/main.cc b/tests_cpp/eigen_2d_euler_sedov_implicit/main.cc
@@ -34,9 +34,8 @@ int main()
   using lin_solver_t = pressio::linearsolvers::Solver<
     pressio::linearsolvers::iterative::Bicgstab, jacob_t>;
   lin_solver_t linSolverObj;
-  auto NonLinSolver=
-    pressio::nonlinearsolvers::create_newton_raphson(stepperObj, linSolverObj);
-  NonLinSolver.setTolerance(1e-5);
+  auto NonLinSolver=pressio::create_newton_solver(stepperObj, linSolverObj);
+  NonLinSolver.setStopTolerance(1e-5);
 
   FomObserver<state_t> Obs("sedov2d_solution.bin", 5);
 

diff --git a/tests_cpp/eigen_2d_euler_sedov_symmetry_implicit/main.cc b/tests_cpp/eigen_2d_euler_sedov_symmetry_implicit/main.cc
@@ -34,9 +34,8 @@ int main()
   using lin_solver_t = pressio::linearsolvers::Solver<
     pressio::linearsolvers::iterative::Bicgstab, jacob_t>;
   lin_solver_t linSolverObj;
-  auto NonLinSolver=
-   pressio::nonlinearsolvers::create_newton_raphson(stepperObj, linSolverObj);
-  NonLinSolver.setTolerance(1e-5);
+  auto NonLinSolver=pressio::create_newton_solver(stepperObj, linSolverObj);
+  NonLinSolver.setStopTolerance(1e-5);
 
   FomObserver<state_t> Obs("sedov2dsym_solution.bin", 5);
 

diff --git a/tests_cpp/eigen_2d_euler_smooth_implicit/main.cc b/tests_cpp/eigen_2d_euler_smooth_implicit/main.cc
@@ -34,9 +34,8 @@ int main()
   using lin_solver_t = pressio::linearsolvers::Solver<
     pressio::linearsolvers::iterative::Bicgstab, jacob_t>;
   lin_solver_t linSolverObj;
-  auto NonLinSolver=
-    pressio::nonlinearsolvers::create_newton_raphson(stepperObj, linSolverObj);
-  NonLinSolver.setTolerance(1e-10);
+  auto NonLinSolver=pressio::create_newton_solver(stepperObj, linSolverObj);
+  NonLinSolver.setStopTolerance(1e-10);
 
   FomObserver<state_t> Obs("smooth2d_solution.bin", 5);
 

diff --git a/tests_cpp/pressio/include/pressio/expressions/as_diagonal_matrix.hpp b/tests_cpp/pressio/include/pressio/expressions/as_diagonal_matrix.hpp
@@ -46,8 +46,8 @@
 //@HEADER
 */
 
-#ifndef EXPRESSIONS_PUBLIC_AS_DIAG_MATRIX_FUNCTIONS_HPP_
-#define EXPRESSIONS_PUBLIC_AS_DIAG_MATRIX_FUNCTIONS_HPP_
+#ifndef EXPRESSIONS_AS_DIAGONAL_MATRIX_HPP_
+#define EXPRESSIONS_AS_DIAGONAL_MATRIX_HPP_
 
 #include "impl/asdiagonalmatrix_traits.hpp"
 #include "impl/asdiagonalmatrix_classes.hpp"
@@ -66,4 +66,4 @@ auto as_diagonal_matrix(T & operand)
 }
 
 }
-#endif  // EXPRESSIONS_PUBLIC_FUNCTIONS_HPP_
+#endif  // EXPRESSIONS_AS_DIAGONAL_MATRIX_HPP_
diff --git a/tests_cpp/pressio/include/pressio/expressions/diag.hpp b/tests_cpp/pressio/include/pressio/expressions/diag.hpp
@@ -46,8 +46,8 @@
 //@HEADER
 */
 
-#ifndef EXPRESSIONS_PUBLIC_DIAG_FUNCTIONS_HPP_
-#define EXPRESSIONS_PUBLIC_DIAG_FUNCTIONS_HPP_
+#ifndef EXPRESSIONS_DIAG_HPP_
+#define EXPRESSIONS_DIAG_HPP_
 
 #include "impl/diag_traits.hpp"
 #include "impl/diag_classes.hpp"
@@ -66,4 +66,4 @@ auto diag(T & operand)
 }
 
 }
-#endif  // EXPRESSIONS_PUBLIC_FUNCTIONS_HPP_
+#endif  // EXPRESSIONS_DIAG_HPP_
diff --git a/tests_cpp/pressio/include/pressio/expressions/impl/asdiagonalmatrix_classes.hpp b/tests_cpp/pressio/include/pressio/expressions/impl/asdiagonalmatrix_classes.hpp
@@ -89,10 +89,8 @@ struct AsDiagonalMatrixExpr<
   {}
 
 public:
-  size_t extent(size_t i) const{
-    assert(i==0 or i==1);
-    (void) i;
-    return extent_;
+  std::size_t extent(size_t i) const{
+    return (i < 2) ? extent_ : std::size_t(1);
   }
 
   native_expr_t const & native() const{

diff --git a/tests_cpp/pressio/include/pressio/expressions/impl/diag_classes.hpp b/tests_cpp/pressio/include/pressio/expressions/impl/diag_classes.hpp
@@ -94,9 +94,7 @@ struct DiagExpr<
   }
 
   size_t extent(size_t i) const{
-    assert(i==0);
-    (void) i;
-    return extent_;
+    return (i < 1) ? extent_ : std::size_t(1);
   }
 
   native_expr_t const & native() const{
@@ -167,9 +165,7 @@ struct DiagExpr<
 
 public:
   size_t extent(size_t i) const{
-    assert(i==0);
-    (void) i;
-    return extent_;
+    return (i < 1) ? extent_ : std::size_t(1);
   }
 
   native_expr_t const & native() const{
@@ -182,7 +178,7 @@ struct DiagExpr<
 
   template<typename _MatrixType = MatrixType>
   mpl::enable_if_t<
-    std::is_same<typename MatrixType::memory_space, Kokkos::HostSpace>::value,
+    std::is_same<typename _MatrixType::memory_space, Kokkos::HostSpace>::value,
     ref_t
     >
   operator()(size_t i) const