Add DG discretizations for LRM, LRMP and GRM

cadet · Apr 23, 2024 · 34694f0 · 34694f0
1 parent 59806c0
commit 34694f0
Show file tree

Hide file tree

Showing 36 changed files with 17,532 additions and 38 deletions.
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
@@ -35,7 +35,8 @@ jobs:
           path: |
             ${{ github.workspace }}/hdf5
             ${{ github.workspace }}/suitesparse/install
-          key: ${{ runner.os }}-deps-${{ steps.get-msvc.outputs.version }}
+            ${{ github.workspace }}/eigen/install
+          key: ${{ runner.os }}-deps-${{ steps.get-msvc.outputs.version }}-1
       - name: Download MKL and TBB
         run: |
           cd "${env:BASE_DIR}"
@@ -44,7 +45,7 @@ jobs:
           nuget install intelmkl.static.win-x64 -Version 2023.0.0.25930
           Invoke-WebRequest -Uri "https://gitlab.com/libeigen/eigen/-/archive/master/eigen-master.zip" -OutFile eigen.zip
           7z x eigen.zip
-      - name: Build UMFPACK and HDF5
+      - name: Build UMFPACK and HDF5 and Eigen3
         if: steps.cache.outputs.cache-hit != 'true'
         run: |
           $base_dir = $($env:BASE_DIR.Replace('\', '/'))
@@ -69,6 +70,14 @@ jobs:
           $ENV:MKLROOT="${env:BASE_DIR}/intelmkl.static.win-x64.2023.0.0.25930/lib/native/win-x64".Replace('\', '/')
           cmake -DCMAKE_INSTALL_PREFIX="${base_dir}\suitesparse\install" -DBLA_VENDOR=Intel10_64lp_seq -DBLA_STATIC=ON -G "Ninja" -DCMAKE_C_FLAGS="/GL" -DCMAKE_STATIC_LINKER_FLAGS="/LTCG" -DCMAKE_BUILD_TYPE=Release -DBUILD_METIS=OFF ..\suitesparse-metis-for-windows-e8d953dffb8a99aa8b65ff3ff03e12a3ed72f90c\
           ninja install
+          cd "${env:BASE_DIR}"
+          Invoke-WebRequest -Uri "https://gitlab.com/libeigen/eigen/-/archive/3.4.0/eigen-3.4.0.zip" -OutFile eigen.zip
+          7z x eigen.zip -oeigen\code -y
+          cd eigen
+          mkdir build
+          cd build
+          cmake ../code/eigen-3.4.0 -DCMAKE_INSTALL_PREFIX="../install"
+          cmake --build . --target install
       - name: Build and Install
         run: |
           cd "${env:BASE_DIR}"
@@ -78,6 +87,7 @@ jobs:
           $ENV:MKLROOT="${env:BASE_DIR}/intelmkl.static.win-x64.2023.0.0.25930/lib/native/win-x64".Replace('\', '/')
           $ENV:TBB_ROOT="${env:BASE_DIR}/inteltbb.devel.win.2021.8.0.25874/lib/native".Replace('\', '/')
           $ENV:UMFPACK_ROOT="${env:BASE_DIR}/suitesparse/install".Replace('\', '/')
+          $ENV:Eigen3_DIR="${env:BASE_DIR}/eigen/install".Replace('\', '/')
           $install_prefix = $($env:INSTALL_PREFIX.Replace('\', '/'))
           $src_dir = $($env:SRC_DIR.Replace('\', '/'))
           $base_dir = $($env:BASE_DIR.Replace('\', '/'))

diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -105,6 +105,9 @@ add_feature_info(ENABLE_DEBUG_THREADING ENABLE_DEBUG_THREADING "Use multi-thread
 option(ENABLE_GRM_2D "Build 2D general rate model" ON)
 add_feature_info(ENABLE_GRM_2D ENABLE_GRM_2D "Build 2D general rate model")
 
+option(ENABLE_DG "Build DG variants of models" ON)
+add_feature_info(ENABLE_DG ENABLE_DG "Build DG variants of models")
+
 option(ENABLE_SUNDIALS_OPENMP "Prefer OpenMP vector implementation of SUNDIALS if available (for large problems)" OFF)
 add_feature_info(ENABLE_SUNDIALS_OPENMP ENABLE_SUNDIALS_OPENMP "Prefer OpenMP vector implementation of SUNDIALS if available (for large problems)")
 
@@ -394,6 +397,21 @@ if (ENABLE_GRM_2D)
 	)
 endif()
 
+set(EIGEN_TARGET "")
+if (ENABLE_DG)
+	find_package(Eigen3 3.4 NO_MODULE)
+
+	# Disable DG if Eigen is not present
+	if (NOT TARGET Eigen3::Eigen)
+		message(STATUS "Disabling DG support because Eigen3 could not be found")
+		set(ENABLE_DG OFF)
+	else()
+		set(EIGEN_TARGET "Eigen3::Eigen")
+		get_target_property(Eigen3_INCLUDE_DIRS Eigen3::Eigen INTERFACE_INCLUDE_DIRECTORIES)
+		include_directories(${Eigen3_INCLUDE_DIRS} "${Eigen3_INCLUDE_DIRS}/..")
+	endif()
+endif()
+
 set(IPO_AVAILABLE OFF)
 if (ENABLE_IPO)
 	include(CheckIPOSupported)
@@ -592,6 +610,14 @@ if (ENABLE_GRM_2D)
 	endif()
 endif()
 
+if (ENABLE_DG)
+	message("Found Eigen3: ${Eigen3_FOUND}")
+	if (TARGET Eigen3::Eigen)
+		message("  Version ${Eigen3_VERSION}")
+		message("  Includes ${Eigen3_INCLUDE_DIRS}")
+	endif()
+endif()
+
 message("Found HDF5: ${HDF5_FOUND}")
 if (HDF5_FOUND)
 	message("  Version ${HDF5_VERSION}")

diff --git a/doc/interface/sensitivities.rst b/doc/interface/sensitivities.rst
@@ -39,7 +39,7 @@ Group /input/sensitivity/param_XXX
 
 ``SENS_NAME``
 
-   Name of the parameter
+   Name of the parameter (Note that ``PAR_RADIUS`` and ``PAR_CORE_RADIUS`` sensitivities are only available for Finite Volume discretization)
 
    ================  ===========================
    **Type:** string  **Length:** :math:`\geq 1`

diff --git a/doc/interface/unit_operations/general_rate_model.rst b/doc/interface/unit_operations/general_rate_model.rst
@@ -454,21 +454,28 @@ For information on model equations, refer to :ref:`general_rate_model_model`.
    ================  ========================  =============================================================================
 
 
-Group /input/model/unit_XXX/discretization - UNIT_TYPE - GENERAL_RATE_MODEL
----------------------------------------------------------------------------
 
+Discretization Methods
+----------------------
+
+CADET has two discretization frameworks available, Finite Volumes (FV) and Discontinuous Galerkin (DG), only one needs to be specified. Both methods approximate the same solution to the same underlying model but can differ regarding computational performance.
+
+Group /input/model/unit_XXX/discretization - UNIT_TYPE - GENERAL_RATE_MODEL
+----------------------------------------------------------------------------------------
+Finite Volumes (Default)
+------------------------
 
 ``NCOL``
 
-   Number of axial column discretization cells
+   Number of axial column discretization points
 
    =============  =========================  =============
    **Type:** int  **Range:** :math:`\geq 1`  **Length:** 1
    =============  =========================  =============
 
 ``NPAR``
 
-   Number of particle (radial) discretization cells for each particle type
+   Number of particle (radial) discretization points for each particle type
 
    =============  =========================  =================================================
    **Type:** int  **Range:** :math:`\geq 1`  **Length:** :math:`1` / :math:`\texttt{NPARTYPE}`
@@ -556,3 +563,95 @@ Group /input/model/unit_XXX/discretization - UNIT_TYPE - GENERAL_RATE_MODEL
 
 For further discretization parameters, see also :ref:`flux_restruction_methods`, and :ref:`non_consistency_solver_parameters`.
 
+Group /input/model/unit_XXX/discretization - UNIT_TYPE - GENERAL_RATE_MODEL_DG
+----------------------------------------------------------------------------------------
+Discontinuous Galerkin
+----------------------
+
+``POLYDEG``
+
+   DG polynomial degree. Optional, defaults to 4. The total number of axial discrete points is given by (``POLYDEG`` + 1 ) * ``NCOL``.
+
+   =============  =========================  =============
+   **Type:** int  **Range:** :math:`\geq 1`  **Length:** 1
+   =============  =========================  =============
+
+``NCOL``
+
+   Number of axial column discretization DG cells\elements. The total number of axial discrete points is given by (``POLYDEG`` + 1 ) * ``NCOL``.
+
+   =============  =========================  =============
+   **Type:** int  **Range:** :math:`\geq 1`  **Length:** 1
+   =============  =========================  =============
+
+``EXACT_INTEGRATION``
+
+   Specifies the DG integration method. Optional, defaults to 0: Choose 1 for exact integration (more accurate but slower), 0 for LGL quadrature (less accurate but faster, typically more performant).
+
+   =============  ===========================  =============
+   **Type:** int  **Range:** :math:`\{0, 1\}`  **Length:** 1
+   =============  ===========================  =============
+
+``NPARTYPE``
+
+   Number of particle types. Optional, inferred from the length of :math:`\texttt{NPAR}` or :math:`\texttt{NBOUND}` if left out.
+
+   =============  =========================  =============
+   **Type:** int  **Range:** :math:`\geq 1`  **Length:** 1
+   =============  =========================  =============
+
+``PARPOLYDEG``
+
+   DG particle (radial) polynomial degree. Optional, defaults to 3. The total number of particle (radial) discrete points is given by (``PARPOLYDEG`` + 1 ) * ``NPARCELL``.
+
+   =============  =========================  =============
+   **Type:** int  **Range:** :math:`\geq 1`  **Length:** 1
+   =============  =========================  =============
+
+``NPARCELL``
+
+   Number of particle (radial) discretization DG cells for each particle type. For the particle discretization, it is usually most performant to fix ``NPARCELL`` = 1 and to increase the polynomial degree for more accuracy.
+
+   =============  =========================  =================================================
+   **Type:** int  **Range:** :math:`\geq 1`  **Length:** :math:`1` / :math:`\texttt{NPARTYPE}`
+   =============  =========================  =================================================
+
+``NBOUND``
+
+   Number of bound states for each component
+
+   =============  =========================  ==================================
+   **Type:** int  **Range:** :math:`\geq 0`  **Length:** :math:`\texttt{NCOMP}`
+   =============  =========================  ==================================
+
+``PAR_GEOM``
+
+   Specifies the particle geometry for all or each particle type. Valid values are :math:`\texttt{SPHERE}`, :math:`\texttt{CYLINDER}`, :math:`\texttt{SLAB}`. Optional, defaults to :math:`\texttt{SPHERE}`.
+
+   ================  =================================================
+   **Type:** string  **Length:** :math:`1` / :math:`\texttt{NPARTYPE}`
+   ================  =================================================
+
+``PAR_DISC_TYPE``
+
+   Specifies the discretization scheme inside the particles for all or each particle type. Valid values are :math:`\texttt{EQUIDISTANT_PAR}`, :math:`\texttt{EQUIVOLUME_PAR}`, and :math:`\texttt{USER_DEFINED_PAR}`.
+
+   ================  =================================================
+   **Type:** string  **Length:** :math:`1` / :math:`\texttt{NPARTYPE}`
+   ================  =================================================
+
+``PAR_DISC_VECTOR``
+
+   Node coordinates for the cell boundaries (ignored if :math:`\texttt{PAR_DISC_TYPE} \neq \texttt{USER_DEFINED_PAR}`). The coordinates are relative and have to include the endpoints :math:`0` and :math:`1`. They are later linearly mapped to the true radial range :math:`[r_{c,j}, r_{p,j}]`. The coordinates for each particle type are appended to one long vector in type-major ordering.
+
+   ================  ========================  ================================================
+   **Type:** double  **Range:** :math:`[0,1]`  **Length:** :math:`\sum_i (\texttt{NPAR}_i + 1)`
+   ================  ========================  ================================================
+
+``USE_ANALYTIC_JACOBIAN``
+
+   Determines whether analytically computed Jacobian matrix (faster) is used (value is 1) instead of Jacobians generated by algorithmic differentiation (slower, value is 0)
+
+   =============  ===========================  =============
+   **Type:** int  **Range:** :math:`\{0, 1\}`  **Length:** 1
+   =============  ===========================  =============
diff --git a/doc/interface/unit_operations/lumped_rate_model_with_pores.rst b/doc/interface/unit_operations/lumped_rate_model_with_pores.rst
@@ -247,13 +247,19 @@ For information on model equations, refer to :ref:`lumped_rate_model_with_pores_
    ================  ========================  =======================================================================
 
 
+Discretization Methods
+----------------------
+
+CADET has two discretization frameworks available, Finite Volumes (FV) and Discontinuous Galerkin (DG), only one needs to be specified. Both methods approximate the same solution to the same underlying model but can differ regarding computational performance.
+
 Group /input/model/unit_XXX/discretization - UNIT_TYPE = LUMPED_RATE_MODEL_WITH_PORES
 -------------------------------------------------------------------------------------
+Finite Volumes (Default)
+------------------------
 
-
 ``NCOL``
 
-   Number of axial column discretization cells
+   Number of axial column discretization points
 
    =============  =========================  =============
    **Type:** int  **Range:** :math:`\geq 1`  **Length:** 1
@@ -308,3 +314,49 @@ Group /input/model/unit_XXX/discretization - UNIT_TYPE = LUMPED_RATE_MODEL_WITH_
    ================  =========================  =============
 
 For further discretization parameters, see also :ref:`flux_restruction_methods`, and :ref:`non_consistency_solver_parameters`.
+
+
+Group /input/model/unit_XXX/discretization - UNIT_TYPE = LUMPED_RATE_MODEL_WITH_PORES_DG
+----------------------------------------------------------------------------------------
+Discontinuous Galerkin
+----------------------
+
+``POLYDEG``
+
+   DG polynomial degree. Optional, defaults to 4. The total number of axial discrete points is given by (``POLYDEG`` + 1 ) * ``NCOL``.
+
+   =============  =========================  =============
+   **Type:** int  **Range:** :math:`\geq 1`  **Length:** 1
+   =============  =========================  =============
+
+``NCOL``
+
+   Number of axial column discretization DG cells\elements. The total number of axial discrete points is given by (``POLYDEG`` + 1 ) * ``NCOL``.
+
+   =============  =========================  =============
+   **Type:** int  **Range:** :math:`\geq 1`  **Length:** 1
+   =============  =========================  =============
+
+``EXACT_INTEGRATION``
+
+   Specifies the DG integration method. Optional, defaults to 0: Choose 1 for exact integration (more accurate but slower), 0 for LGL quadrature (less accurate but faster, typically more performant).
+
+   =============  ===========================  =============
+   **Type:** int  **Range:** :math:`\{0, 1\}`  **Length:** 1
+   =============  ===========================  =============
+
+``NBOUND``
+
+   Number of bound states for each component
+
+   =============  =========================  ==================================
+   **Type:** int  **Range:** :math:`\geq 0`  **Length:** :math:`\texttt{NCOMP}`
+   =============  =========================  ==================================
+
+``USE_ANALYTIC_JACOBIAN``
+
+   Determines whether analytically computed Jacobian matrix (faster) is used (value is 1) instead of Jacobians generated by algorithmic differentiation (slower, value is 0)
+
+   =============  ===========================  =============
+   **Type:** int  **Range:** :math:`\{0, 1\}`  **Length:** 1
+   =============  ===========================  =============
diff --git a/doc/interface/unit_operations/lumped_rate_model_without_pores.rst b/doc/interface/unit_operations/lumped_rate_model_without_pores.rst
@@ -141,13 +141,19 @@ For information on model equations, refer to :ref:`lumped_rate_model_without_por
    ================  =====================  =============
 
 
+Discretization Methods
+----------------------
+
+CADET has two discretization frameworks available, Finite Volumes (FV) and Discontinuous Galerkin (DG), only one needs to be specified. Both methods approximate the same solution to the same underlying model but can differ regarding computational performance.
+
 Group /input/model/unit_XXX/discretization - UNIT_TYPE = LUMPED_RATE_MODEL_WITHOUT_PORES
 ----------------------------------------------------------------------------------------
-
+Finite Volumes (Default)
+------------------------
 
 ``NCOL``
 
-   Number of axial column discretization cells
+   Number of axial column discretization points
 
    =============  =========================  =============
    **Type:** int  **Range:** :math:`\geq 1`  **Length:** 1
@@ -163,11 +169,47 @@ Group /input/model/unit_XXX/discretization - UNIT_TYPE = LUMPED_RATE_MODEL_WITHO
 
 ``RECONSTRUCTION``
 
-   Type of reconstruction method for fluxes
+   Type of reconstruction method for fluxes only (only needs to be specified for FV)
 
    ================  ================================  =============
    **Type:** string  **Range:** :math:`\texttt{WENO}`  **Length:** 1
    ================  ================================  =============
 
 For further discretization parameters, see also :ref:`flux_restruction_methods`, and :ref:`non_consistency_solver_parameters`.
 
+Group /input/model/unit_XXX/discretization - UNIT_TYPE = LUMPED_RATE_MODEL_WITHOUT_PORES_DG
+-------------------------------------------------------------------------------------------
+Discontinuous Galerkin
+----------------------
+
+``POLYDEG``
+
+   DG polynomial degree. Optional, defaults to 4. The total number of axial discrete points is given by (``POLYDEG`` + 1 ) * ``NCOL``.
+
+   =============  =========================  =============
+   **Type:** int  **Range:** :math:`\geq 1`  **Length:** 1
+   =============  =========================  =============
+
+``NCOL``
+
+   Number of axial column discretization DG cells\elements. The total number of axial discrete points is given by (``POLYDEG`` + 1 ) * ``NCOL``.
+
+   =============  =========================  =============
+   **Type:** int  **Range:** :math:`\geq 1`  **Length:** 1
+   =============  =========================  =============
+
+``EXACT_INTEGRATION``
+
+   Specifies the DG integration method. Optional, defaults to 0: Choose 1 for exact integration (more accurate but slower), 0 for LGL quadrature (less accurate but faster, typically more performant).
+
+   =============  ===========================  =============
+   **Type:** int  **Range:** :math:`\{0, 1\}`  **Length:** 1
+   =============  ===========================  =============
+
+``USE_ANALYTIC_JACOBIAN``
+
+   Determines whether analytically computed Jacobian matrix (faster) is used (value is 1) instead of Jacobians generated by algorithmic differentiation (slower, value is 0)
+
+   =============  ===========================  =============
+   **Type:** int  **Range:** :math:`\{0, 1\}`  **Length:** 1
+   =============  ===========================  =============
diff --git a/doc/literature.bib b/doc/literature.bib
@@ -436,4 +436,27 @@ @article{Jaepel2022
 url = {https://www.sciencedirect.com/science/article/pii/S0021967322005830},
 author = {Ronald Colin Jäpel and Johannes Felix Buyel},
 }
-
+@book{Kopriva2009,
+address = {Dordrecht},
+author = {Kopriva, David A.},
+series = {Scientific {Computation}},
+title = {Implementing {Spectral} {Methods} for {Partial} {Differential} {Equations}: {Algorithms} for {Scientists} and {Engineers}},
+isbn = {978-90-481-2260-8 978-90-481-2261-5},
+shorttitle = {Implementing {Spectral} {Methods} for {Partial} {Differential} {Equations}},
+url = {http://link.springer.com/10.1007/978-90-481-2261-5},
+urldate = {2024-01-12},
+publisher = {Springer Netherlands},
+year = {2009},
+doi = {https://doi.org/10.1007/978-90-481-2261-5},
+}
+@article{Breuer2023,
+title = {Spatial discontinuous Galerkin spectral element method for a family of chromatography models in CADET},
+journal = {Computers \& Chemical Engineering},
+volume = {177},
+pages = {108340},
+year = {2023},
+issn = {0098-1354},
+doi = {https://doi.org/10.1016/j.compchemeng.2023.108340},
+url = {https://www.sciencedirect.com/science/article/pii/S0098135423002107},
+author = {Jan Michael Breuer and Samuel Leweke and Johannes Schmölder and Gregor Gassner and Eric {von Lieres}},
+}