From a0fcd73531b3f33410f5b8e462b5a16e1c9752d6 Mon Sep 17 00:00:00 2001 From: Jan Breuer <74359367+jbreue16@users.noreply.github.com> Date: Mon, 18 Dec 2023 17:41:11 +0100 Subject: [PATCH] Streamline DG Discretization Interface (#152) * Change DG unit interfaces from unit_type suffix to discretization specification * Infer DG disc from NCOL and NPAR and set default values * Rename DG element fields to nelem and par_nelem * Fix test interfaces, i.e. adjust to createUnitOperation function --- .../spatial_discretization_methods.rst | 125 ++++++++++++++ .../unit_operations/general_rate_model.rst | 153 +++++++---------- .../lumped_rate_model_with_pores.rst | 78 ++++----- .../lumped_rate_model_without_pores.rst | 79 ++++----- include/cadet/ModelBuilder.hpp | 22 +-- src/libcadet/CMakeLists.txt | 3 + src/libcadet/ModelBuilderImpl.cpp | 58 +++---- src/libcadet/ModelBuilderImpl.hpp | 4 +- src/libcadet/model/GeneralRateModel.cpp | 14 +- src/libcadet/model/GeneralRateModel.hpp | 3 + src/libcadet/model/GeneralRateModel2D.cpp | 6 +- .../model/GeneralRateModelBuilder.cpp | 91 ++++++++++ src/libcadet/model/GeneralRateModelDG.cpp | 157 +++++++++--------- src/libcadet/model/GeneralRateModelDG.hpp | 2 +- src/libcadet/model/GeneralRateModelDGFV.cpp | 16 -- src/libcadet/model/GeneralRateModelDGFV.hpp | 2 +- src/libcadet/model/InletModel.cpp | 4 +- .../model/LumpedRateModelWithPores.cpp | 14 +- .../model/LumpedRateModelWithPores.hpp | 3 + .../model/LumpedRateModelWithPoresBuilder.cpp | 87 ++++++++++ .../model/LumpedRateModelWithPoresDG.cpp | 44 +++-- .../model/LumpedRateModelWithoutPores.cpp | 18 +- .../model/LumpedRateModelWithoutPores.hpp | 3 + .../LumpedRateModelWithoutPoresBuilder.cpp | 88 ++++++++++ .../model/LumpedRateModelWithoutPoresDG.cpp | 37 +++-- src/libcadet/model/OutletModel.cpp | 4 +- src/libcadet/model/StirredTankModel.cpp | 4 +- test/CSTR-Residual.cpp | 2 +- test/ColumnTests.cpp | 2 +- test/GeneralRateModel2D.cpp | 32 ++-- test/UnitOperationTests.cpp | 2 +- 31 files changed, 765 insertions(+), 392 deletions(-) create mode 100644 doc/interface/spatial_discretization_methods.rst create mode 100644 src/libcadet/model/GeneralRateModelBuilder.cpp create mode 100644 src/libcadet/model/LumpedRateModelWithPoresBuilder.cpp create mode 100644 src/libcadet/model/LumpedRateModelWithoutPoresBuilder.cpp diff --git a/doc/interface/spatial_discretization_methods.rst b/doc/interface/spatial_discretization_methods.rst new file mode 100644 index 000000000..ed6060520 --- /dev/null +++ b/doc/interface/spatial_discretization_methods.rst @@ -0,0 +1,125 @@ +.. spatial_discretization_methods: + +Spatial discretization methods +============================== + +Group /input/model/unit_XXX/discretization/SPATIAL_METHOD - Details on the methods +---------------------------------------------------------------------------------- + +CADET offers two spatial discretization methods: Finite Volumes (FV) and Discontinuous Galerkin (DG). +While both methods approximate the same solution to the underlying models, they may differ in terms of computational performance. +Generally, FV can be more performant for small problem sizes and solutions with steep gradients, while DG excels for large problem sizes and smooth solutions. + +In the following, we give a brief introduction to the numerical theory that is most relevant for the computational performance of the methods. +Based on that theory and our experience, we give advice on which method to use in which scenario, how to identify the more performant method, and how to specify the discretization parameters. +For a comprehensive description on the FV and DG methods as they are implemented in CADET, we refer to our publications on `CADET-FV `_ and `CADET-DG `_. + +Discrete system size +-------------------- + +Numerical methods discretize the continuous (here: spatial) domain of the equations into a finite set of discrete points. +Then, a system of equations is formulated for those points, and both system size and number of unknowns / degrees of freedom (DoF) are given by the number of discrete points. +This system can be linear or non-linear, depending on the method (this will become important again in the section on smooth solutions). +The wall clock time to compute the solution depends on the system size and, thus, on the number of discrete points. +Conversely, the numerical solution is more accurate with more discrete points. +Thus, we trade computation time for approximation accuracy by specifying the parameters that determine the number of discrete points. +For the FV scheme, the number of axial discrete points in the column is given by the number of volume cells ``NCOL``. +For the DG scheme, the number of axial discrete points in the column is given by the number of polynomial interpolation nodes (= ``POLYDEG`` + 1) times the number of DG elements ``NELEM``. +The LRMP and GRM additionally consider particle equations that are also discretized. +In the spatially discretized equations, a single particle is incorporated at each axial discrete point, which increases the total number of DOF per axial point, especially for the GRM where particles are fully resolved. +The parameters for the GRM particle discretization are given for FV in ``NPAR`` and for DG in ``PAR_POLYDEG`` and ``PAR_NELEM``. + +Order of convergence +-------------------- + +The computational performance of a numerical method depends on its theoretical order of convergence. +The order of convergence refers to the rate at which the method's approximation approaches the exact solution under refinement of the spatial grid. + +A higher-order method can be faster than a low-order method: +Imagine a high- and a low-order method's approximation to exhibit similarly bad approximation accuracy due to a coarse spatial resolution. +Refining the grid for both methods by the same number of discrete points improves the approximation accuracy of the higher-order method more than the other one. +Thus, the low-order method requires more DOFs and ultimately more compute time to compute a solution of the same accuracy. + +The theoretical order of convergence is an asymptotic property, however. +Having the exact solution, we can compute an experimental order of convergence (EOC) via the formula + +.. math:: + :label: ModelColumn + + \begin{aligned} + EOC_k = \frac{log(\varepsilon_{k+1} / \varepsilon_{k}) }{log(n_{k} / n_{k+1})}, + \end{aligned} + +with :math:`\varepsilon_{k}` and :math:`n_{k}` denoting some error norm and the degrees of freedom of the kth approximation. +The EOC approaches the theoretical order of convergence for :math:`k \rightarrow \infty` but is typically lower for underresolved problems. +High-order methods typically suffer from start-off problems, i.e. they typically won't exhibit their high order until the grid is fine enough and a certain accuracy is already reached. +That is, increasing the number of discrete points from, e.g., 2 to 4 typically does not improve the solution according to the theoretical order of convergence but by a much smaller EOC. +The EOC is highly problem-dependent, and it is generally unknown when a high-order method will actually be faster than a lower-order method. +Experience shows that higher-order methods work well for smooth solutions. + +The theoretical order of convergence for the CADET-FV scheme is fixed at 2. +For the CADET-DG scheme, it is :math:`N_d + 1` with :math:`N_d` denoting the polynomial degree, and can thus be user-defined by specifying the field ``POLYDEG`` (and ``PAR_POLYDEG`` for the GRM). +As a convergence order of :math:`\gt 6` is hardly realized within the approximation error of engineering tolerance (due to start-off problems), we recommend a maximum polynomial order of 5. +As the FV scheme oftentimes yields an EOC of around 2.5 and is computationally more enhanced (less arithmetic operations per DOF and customized factorization) than the DG code, we recommend a polynomial degree of at least 3 to top this. + +Smooth solutions +---------------- + +Smoothness in functions is characterized by the absence of sudden changes, reflected in the continuity and differentiability of the function and its derivatives. +In numerical simulation, smoothness indicators are often based on derivatives or moments of the solution. +That is, strong gradients and high frequencies are used to identify non-smooth parts of the solution. +Godunov's order barrier theorem shows why the concept of smoothness plays a crucial role in the deployment of numerical methods. +It states that linear high-order methods that are monotonous are at most first-order accurate. +Linear higher-order (:math:`\gt 1`) methods thus suffer from artificial oscillations at non-smooth parts of the solution, specifically at discontinuities and strong gradients. +Some higher-order methods, such as CADET-FV (2nd order), contain a non-linear mechanism to suppress these oscillations. +The non-linear WENO mechanism employed in CADET-FV can be fine-tuned via the fields specified here :ref:`flux_restruction_methods`. +Unfortunately, non-linear higher-order methods (order :math:`\geq 3`) are either not applicable (e.g., undefined boundary treatment) or have other shortcomings, such as more highly problem-dependent parameters. + +CADET-DG is a linear high-order method (arbitrary order) and thus exhibits oscillatory behaviour at strong gradients, which increases the approximation error and results in a smaller EOC for lower resolutions. +As strong gradients are a local phenomena that can be captured by employing more discrete points, DG becomes more performant again for higher resolutions. +This, however, might happen after the engineering error tolerance is by far surpassed. +Hence, CADET-FV as a stabilized lower-order method can be more performant, depending on the setting. +The DG scheme reduces its oscillatory behaviour by adding artificial numerical dispersion at element interfaces. +Thus, the use of a lower polynomial degree and more elements is recommended for rather non-smooth problems. + +In Chromatography, mathematical discontinuities never happen, as there are always some dispersive effects in reality. +Chromatography models, however, allow for discontinuities if dispersion parameters are set to zero. +Moreover, steep and self-sharpening concentration fronts might appear due to binding. +Binding models that might cause self-sharpening concentration fronts are often associated with competitive Langmuir type isotherms for components with differently strong binding properties. +Nonetheless, a lot of chromatography settings yield rather smooth concentration profiles, for which DG is the better choice in terms of computational performance. + +Recommendations on the choice of spatial discretization methods +--------------------------------------------------------------- + +We recommend the FV method for + +- Small problem sizes, e.g., low spatial resolution with the LRM +- Problems with strong gradients, e.g., no or low dispersion and bindings that create sharp fronts +- Bindings that mathematically require positive values or exhibit strange behaviour with negative concentration values + +We recommend the DG method for + +- Large problem sizes, e.g., high resolutions and more complex models (i.e. the LRMP and specifically the GRM) +- Smooth problems, e.g., sufficient dispersion + +Recommendations on DG discretization parameters +----------------------------------------------- + +- Employ an axial polynomial degree between 3 and 5 +- Select a lower axial polynomial degree for non smooth tendency and employ more elements instead. Converse choice for smooth problems +- Adjust the DG particle polynomial degree to control approximation accuracy; leave the number of elements at one. Make exceptions if very steep gradients occur inside the particles or when specific parts of the particle domain are more interesting (here, you can resolve more interesting parts by a user-defined spacing of multiple elements) +- The field ``EXACT_INTEGRATION`` specifies the DG polynomial integration method. The default value of 0 (collocation DG) is expected to be slightly more performant in most settings + +Refinement strategy +------------------- + +A common problem in numerical simulation is that the number of discrete points required to yield an accurate approximation within a specific tolerance is unknown. +We thus recommend determining the approximation error via comparison with a refined reference approximation. +Both the theoretical order of convergence and the EOC can be used to estimate the required number of discrete points. + +Note on DG solution vector +-------------------------- + +Any liquid or solid concentration within the column or particles is reported on the discrete points that are employed by the method. +That is, DG yields a piece-wise polynomial approximation on Lagrange-Gauss-Lobatto nodes. +If the solution is desired on a different grid, element-wise polynomial interpolation should be applied, and element interface values must be averaged. diff --git a/doc/interface/unit_operations/general_rate_model.rst b/doc/interface/unit_operations/general_rate_model.rst index 6528131d0..5f55adef2 100644 --- a/doc/interface/unit_operations/general_rate_model.rst +++ b/doc/interface/unit_operations/general_rate_model.rst @@ -426,26 +426,6 @@ Group /input/model/unit_XXX - UNIT_TYPE - GENERAL_RATE_MODEL **Type:** double **Range:** :math:`[0,1]` **Length:** :math:`\texttt{NPARTYPE} / \texttt{NCOL} \cdot \texttt{NPARTYPE}` ================ ======================== ============================================================================= - - -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 points - - ============= ========================= ============= - **Type:** int **Range:** :math:`\geq 1` **Length:** 1 - ============= ========================= ============= - ``NPARTYPE`` Number of particle types. Optional, inferred from the length of :math:`\texttt{NPAR}` or :math:`\texttt{NBOUND}` if left out. @@ -454,14 +434,6 @@ Finite Volumes (Default) **Type:** int **Range:** :math:`\geq 1` **Length:** 1 ============= ========================= ============= -``NPAR`` - - Number of particle (radial) discretization points for each particle type - - ============= ========================= ================================================= - **Type:** int **Range:** :math:`\geq 1` **Length:** :math:`1` / :math:`\texttt{NPARTYPE}` - ============= ========================= ================================================= - ``NBOUND`` Number of bound states for each component in each particle type in type-major ordering @@ -494,14 +466,6 @@ Finite Volumes (Default) **Type:** double **Range:** :math:`[0,1]` **Length:** :math:`\sum_i (\texttt{NPAR}_i + 1)` ================ ======================== ================================================ -``PAR_BOUNDARY_ORDER`` - - Order of accuracy of outer particle boundary condition. Optional, defaults to :math:`2`. - - ============= ============================ ============= - **Type:** int **Range:** :math:`\{ 1,2 \}` **Length:** 1 - ============= ============================ ============= - ``USE_ANALYTIC_JACOBIAN`` Determines whether analytically computed Jacobian matrix (faster) is used (value is :math:`1`) instead of Jacobians generated by algorithmic differentiation (slower, value is :math:`0`) @@ -510,6 +474,49 @@ Finite Volumes (Default) **Type:** int **Range:** :math:`\{0, 1\}` **Length:** 1 ============= =========================== ============= +Spatial discretization - Numerical Methods +------------------------------------------ + +CADET offers two spatial discretization methods: Finite Volumes (FV) and Discontinuous Galerkin (DG). Only one method needs to be specified. +While both methods approximate the same solution to the same underlying model, they may differ in terms of computational performance. +Generally, FV is more performant for solutions with steep gradients, while DG excels for smooth solutions. +For further information on the choice of discretization methods and their parameters, see :ref:`spatial_discretization_methods`. + +``SPATIAL_METHOD`` + + Spatial discretization method. Optional, defaults to :math:`\texttt{FV}` + + ================ =============================================== ============= + **Type:** string **Range:** :math:`\{\texttt{FV}, \texttt{DG}\}` **Length:** 1 + ================ =============================================== ============= + +Finite Volumes (Default) +------------------------ + +``NCOL`` + + Number of axial column discretization points + + ============= ========================= ============= + **Type:** int **Range:** :math:`\geq 1` **Length:** 1 + ============= ========================= ============= + +``NPAR`` + + Number of particle (radial) discretization points for each particle type + + ============= ========================= ================================================= + **Type:** int **Range:** :math:`\geq 1` **Length:** :math:`1` / :math:`\texttt{NPARTYPE}` + ============= ========================= ================================================= + +``PAR_BOUNDARY_ORDER`` + + Order of accuracy of outer particle boundary condition. Optional, defaults to :math:`2`. + + ============= ============================ ============= + **Type:** int **Range:** :math:`\{ 1,2 \}` **Length:** 1 + ============= ============================ ============= + ``RECONSTRUCTION`` Type of reconstruction method for fluxes @@ -558,16 +565,22 @@ Finite Volumes (Default) **Type:** int **Range:** :math:`\{0, 1\}` **Length:** 1 ============= =========================== ============= -For further discretization parameters, see also :ref:`flux_restruction_methods`, and :ref:`non_consistency_solver_parameters`. +For further FV specific discretization parameters, see also :ref:`flux_restruction_methods`. -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``. + DG polynomial degree. Optional, defaults to 4 and :math:`N_d \in \{3, 4, 5\}` is recommended. The total number of axial discrete points is given by (``POLYDEG`` + 1 ) * ``NELEM`` + + ============= ========================= ============= + **Type:** int **Range:** :math:`\geq 1` **Length:** 1 + ============= ========================= ============= + +``NELEM`` + + Number of axial column discretization DG cells\elements. The total number of axial discrete points is given by (``POLYDEG`` + 1 ) * ``NELEM`` ============= ========================= ============= **Type:** int **Range:** :math:`\geq 1` **Length:** 1 @@ -575,7 +588,7 @@ Discontinuous Galerkin ``NCOL`` - Number of axial column discretization DG cells\elements. The total number of axial discrete points is given by (``POLYDEG`` + 1 ) * ``NCOL``. + Number of axial discrete points. Will be ignored if ``NELEM`` is defined, otherwise number of elements is calculated via ``NELEM`` = :math:`\lfloor ``NCOL`` / (``POLYDEG`` + 1 ) \rfloor` and ``NCOL`` is internally overwritten accordingly ============= ========================= ============= **Type:** int **Range:** :math:`\geq 1` **Length:** 1 @@ -583,72 +596,24 @@ Discontinuous Galerkin ``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). + Specifies the DG integration variant. Optional, defaults to 0 ============= =========================== ============= **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`` +``PAR_POLYDEG`` - DG particle (radial) polynomial degree. Optional, defaults to 3. The total number of particle (radial) discrete points is given by (``PARPOLYDEG`` + 1 ) * ``NPARCELL``. + DG particle (radial) polynomial degree. Optional, defaults to 3. The total number of particle (radial) discrete points is given by (``PARPOLYDEG`` + 1 ) * ``PAR_NELEM``. ============= ========================= ============= **Type:** int **Range:** :math:`\geq 1` **Length:** 1 ============= ========================= ============= -``NPARCELL`` +``PAR_NELEM`` - 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. + Number of particle (radial) discretization DG elements for each particle type. For the particle discretization, it is usually most performant to fix ``PAR_NELEM`` = 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 index 9319427f2..209b92529 100644 --- a/doc/interface/unit_operations/lumped_rate_model_with_pores.rst +++ b/doc/interface/unit_operations/lumped_rate_model_with_pores.rst @@ -220,25 +220,9 @@ Group /input/model/unit_XXX - UNIT_TYPE = LUMPED_RATE_MODEL_WITH_PORES **Type:** double **Range:** :math:`[0,1]` **Length:** :math:`\texttt{NPARTYPE}` / :math:`\texttt{NCOL} \cdot \texttt{NPARTYPE}` ================ ======================== ======================================================================= - -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 points - - ============= ========================= ============= - **Type:** int **Range:** :math:`\geq 1` **Length:** 1 - ============= ========================= ============= - ``NPARTYPE`` Number of particle types. Optional, inferred from the length of :math:`\texttt{NBOUND}` if left out. @@ -270,6 +254,33 @@ Finite Volumes (Default) ============= =========================== ============= **Type:** int **Range:** :math:`\{0, 1\}` **Length:** 1 ============= =========================== ============= + +Spatial discretization - Numerical Methods +------------------------------------------ + +``SPATIAL_METHOD`` + + Spatial discretization method. Optional, defaults to :math:`\texttt{FV}` + + ================ =============================================== ============= + **Type:** string **Range:** :math:`\{\texttt{FV}, \texttt{DG}\}` **Length:** 1 + ================ =============================================== ============= + +CADET offers two spatial discretization methods: Finite Volumes (FV) and Discontinuous Galerkin (DG). Only one method needs to be specified. +While both methods approximate the same solution to the same underlying model, they may differ in terms of computational performance. +Generally, FV is more performant for solutions with steep gradients, while DG excels for smooth solutions. +For further information on the choice of discretization methods and their parameters, see :ref:`spatial_discretization_methods`. + +Finite Volumes (Default) +------------------------ + +``NCOL`` + + Number of axial column discretization points + + ============= ========================= ============= + **Type:** int **Range:** :math:`\geq 1` **Length:** 1 + ============= ========================= ============= ``RECONSTRUCTION`` @@ -311,50 +322,41 @@ Finite Volumes (Default) **Type:** double **Range:** :math:`\geq 0` **Length:** 1 ================ ========================= ============= -For further discretization parameters, see also :ref:`flux_restruction_methods`, and :ref:`non_consistency_solver_parameters`. +For further FV specific discretization parameters, see also :ref:`flux_restruction_methods`. - -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``. + DG polynomial degree. Optional, defaults to 4 and :math:`N_d \in \{3, 4, 5\}` is recommended. The total number of axial discrete points is given by (``POLYDEG`` + 1 ) * ``NELEM`` ============= ========================= ============= **Type:** int **Range:** :math:`\geq 1` **Length:** 1 ============= ========================= ============= -``NCOL`` +``NELEM`` - Number of axial column discretization DG cells\elements. The total number of axial discrete points is given by (``POLYDEG`` + 1 ) * ``NCOL``. + Number of axial column discretization DG cells\elements. The total number of axial discrete points is given by (``POLYDEG`` + 1 ) * ``NELEM`` ============= ========================= ============= **Type:** int **Range:** :math:`\geq 1` **Length:** 1 ============= ========================= ============= -``EXACT_INTEGRATION`` +``NCOL`` - 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). + Number of axial discrete points. Will be ignored if ``NELEM`` is defined, otherwise number of elements is calculated via ``NELEM`` = :math:`\lfloor ``NCOL`` / (``POLYDEG`` + 1 ) \rfloor` and ``NCOL`` is internally overwritten accordingly - ============= =========================== ============= - **Type:** int **Range:** :math:`\{0, 1\}` **Length:** 1 - ============= =========================== ============= - -``NBOUND`` + ============= ========================= ============= + **Type:** int **Range:** :math:`\geq 1` **Length:** 1 + ============= ========================= ============= - Number of bound states for each component - - ============= ========================= ================================== - **Type:** int **Range:** :math:`\geq 0` **Length:** :math:`\texttt{NCOMP}` - ============= ========================= ================================== - -``USE_ANALYTIC_JACOBIAN`` +``EXACT_INTEGRATION`` - Determines whether analytically computed Jacobian matrix (faster) is used (value is 1) instead of Jacobians generated by algorithmic differentiation (slower, value is 0) + Specifies the DG integration variant. Optional, defaults to 0 ============= =========================== ============= **Type:** int **Range:** :math:`\{0, 1\}` **Length:** 1 ============= =========================== ============= + +For further general discretization parameters, see also :ref:`non_consistency_solver_parameters`. diff --git a/doc/interface/unit_operations/lumped_rate_model_without_pores.rst b/doc/interface/unit_operations/lumped_rate_model_without_pores.rst index 0755b7e3d..b27030469 100644 --- a/doc/interface/unit_operations/lumped_rate_model_without_pores.rst +++ b/doc/interface/unit_operations/lumped_rate_model_without_pores.rst @@ -128,25 +128,9 @@ Group /input/model/unit_XXX - UNIT_TYPE = LUMPED_RATE_MODEL_WITHOUT_PORES ================ ===================== ============= **Type:** double **Range:** :math:`>0` **Length:** 1 ================ ===================== ============= - - -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 points - - ============= ========================= ============= - **Type:** int **Range:** :math:`\geq 1` **Length:** 1 - ============= ========================= ============= ``NBOUND`` @@ -164,57 +148,76 @@ Finite Volumes (Default) **Type:** int **Range:** :math:`\{0, 1\}` **Length:** 1 ============= =========================== ============= +Spatial discretization - Numerical Methods +------------------------------------------ + +CADET offers two spatial discretization methods: Finite Volumes (FV) and Discontinuous Galerkin (DG). Only one method needs to be specified. +While both methods approximate the same solution to the same underlying model, they may differ in terms of computational performance. +Generally, FV is more performant for solutions with steep gradients, while DG excels for smooth solutions. +For further information on the choice of discretization methods and their parameters, see :ref:`spatial_discretization_methods`. + +``SPATIAL_METHOD`` + + Spatial discretization method. Optional, defaults to :math:`\texttt{FV}` + + ================ =============================================== ============= + **Type:** string **Range:** :math:`\{\texttt{FV}, \texttt{DG}\}` **Length:** 1 + ================ =============================================== ============= + +Finite Volumes (Default) +------------------------ + +``NCOL`` + + Number of axial column discretization points + + ============= ========================= ============= + **Type:** int **Range:** :math:`\geq 1` **Length:** 1 + ============= ========================= ============= + ``RECONSTRUCTION`` - Type of reconstruction method for fluxes only (only needs to be specified for FV) + Type of reconstruction method for fluxes only ================ ================================ ============= **Type:** string **Range:** :math:`\texttt{WENO}` **Length:** 1 ================ ================================ ============= -For further discretization parameters, see also :ref:`flux_restruction_methods`, and :ref:`non_consistency_solver_parameters`. +For further FV specific discretization parameters, see also :ref:`flux_restruction_methods`. -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``. + DG polynomial degree. Optional, defaults to 4 and :math:`N_d \in \{3, 4, 5\}` is recommended. The total number of axial discrete points is given by (``POLYDEG`` + 1 ) * ``NELEM`` ============= ========================= ============= **Type:** int **Range:** :math:`\geq 1` **Length:** 1 ============= ========================= ============= -``NCOL`` +``NELEM`` - Number of axial column discretization DG cells\elements. The total number of axial discrete points is given by (``POLYDEG`` + 1 ) * ``NCOL``. + Number of axial column discretization DG cells\elements. The total number of axial discrete points is given by (``POLYDEG`` + 1 ) * ``NELEM`` ============= ========================= ============= **Type:** int **Range:** :math:`\geq 1` **Length:** 1 ============= ========================= ============= -``EXACT_INTEGRATION`` +``NCOL`` - 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). + Number of axial discrete points. Will be ignored if ``NELEM`` is defined, otherwise number of elements is calculated via ``NELEM`` = :math:`\lfloor ``NCOL`` / (``POLYDEG`` + 1 ) \rfloor` and ``NCOL`` is internally overwritten accordingly - ============= =========================== ============= - **Type:** int **Range:** :math:`\{0, 1\}` **Length:** 1 - ============= =========================== ============= - -``NBOUND`` + ============= ========================= ============= + **Type:** int **Range:** :math:`\geq 1` **Length:** 1 + ============= ========================= ============= - Number of bound states for each component - - ============= ========================= ================================== - **Type:** int **Range:** :math:`\geq 0` **Length:** :math:`\texttt{NCOMP}` - ============= ========================= ================================== - -``USE_ANALYTIC_JACOBIAN`` +``EXACT_INTEGRATION`` - Determines whether analytically computed Jacobian matrix (faster) is used (value is 1) instead of Jacobians generated by algorithmic differentiation (slower, value is 0) + Specifies the DG integration variant. Optional, defaults to 0 ============= =========================== ============= **Type:** int **Range:** :math:`\{0, 1\}` **Length:** 1 ============= =========================== ============= + +For further general discretization parameters, see also :ref:`non_consistency_solver_parameters`. diff --git a/include/cadet/ModelBuilder.hpp b/include/cadet/ModelBuilder.hpp index d05b57617..dfb0d42b7 100644 --- a/include/cadet/ModelBuilder.hpp +++ b/include/cadet/ModelBuilder.hpp @@ -118,17 +118,17 @@ class CADET_API IModelBuilder */ virtual IModel* createUnitOperation(IParameterProvider& paramProvider, UnitOpIdx uoId) = 0; - /** - * @brief Creates a unit operation model - * @details The created unit operation model is not owned by the IModelBuilder. - * Ownership of unit operation models is handled by IModelSystem. - * Unit operation models can be destroyed by calling destroyUnitOperation(). - * - * @param [in] uoType Name of the unit operation model - * @param [in] uoId Unit operation index assigned to the created unit operation - * @return Uninitialized unit operation model or @c nullptr if an error occurred - */ - virtual IModel* createUnitOperation(const std::string& uoType, UnitOpIdx uoId) = 0; + ///** + // * @brief Creates a unit operation model + // * @details The created unit operation model is not owned by the IModelBuilder. + // * Ownership of unit operation models is handled by IModelSystem. + // * Unit operation models can be destroyed by calling destroyUnitOperation(). + // * + // * @param [in] uoType Name of the unit operation model + // * @param [in] uoId Unit operation index assigned to the created unit operation + // * @return Uninitialized unit operation model or @c nullptr if an error occurred + // */ + //virtual IModel* createUnitOperation(const std::string& uoType, UnitOpIdx uoId) = 0; /** * @brief Destroys the given IModel diff --git a/src/libcadet/CMakeLists.txt b/src/libcadet/CMakeLists.txt index ef955287c..fb2ca4623 100644 --- a/src/libcadet/CMakeLists.txt +++ b/src/libcadet/CMakeLists.txt @@ -132,8 +132,11 @@ set(LIBCADET_SOURCES ${CMAKE_SOURCE_DIR}/src/libcadet/model/OutletModel.cpp ${CMAKE_SOURCE_DIR}/src/libcadet/model/StirredTankModel.cpp ${CMAKE_SOURCE_DIR}/src/libcadet/model/LumpedRateModelWithoutPores.cpp + ${CMAKE_SOURCE_DIR}/src/libcadet/model/LumpedRateModelWithoutPoresBuilder.cpp ${CMAKE_SOURCE_DIR}/src/libcadet/model/LumpedRateModelWithPores.cpp + ${CMAKE_SOURCE_DIR}/src/libcadet/model/LumpedRateModelWithPoresBuilder.cpp ${CMAKE_SOURCE_DIR}/src/libcadet/model/GeneralRateModel.cpp + ${CMAKE_SOURCE_DIR}/src/libcadet/model/GeneralRateModelBuilder.cpp ${CMAKE_SOURCE_DIR}/src/libcadet/model/ParameterMultiplexing.cpp ${CMAKE_SOURCE_DIR}/src/libcadet/model/parts/ConvectionDispersionOperator.cpp ${CMAKE_SOURCE_DIR}/src/libcadet/model/parts/ConvectionDispersionOperatorDG.cpp diff --git a/src/libcadet/ModelBuilderImpl.cpp b/src/libcadet/ModelBuilderImpl.cpp index 6f420257b..5bcd33d19 100644 --- a/src/libcadet/ModelBuilderImpl.cpp +++ b/src/libcadet/ModelBuilderImpl.cpp @@ -30,19 +30,15 @@ namespace cadet { namespace model { - void registerInletModel(std::unordered_map>& models); - void registerOutletModel(std::unordered_map>& models); - - void registerGeneralRateModel(std::unordered_map>& models); - void registerGeneralRateModelDG(std::unordered_map>& models); - void registerGeneralRateModelDGFV(std::unordered_map>& models); - void registerLumpedRateModelWithPores(std::unordered_map>& models); - void registerLumpedRateModelWithPoresDG(std::unordered_map>& models); - void registerLumpedRateModelWithoutPores(std::unordered_map>& models); - void registerLumpedRateModelWithoutPoresDG(std::unordered_map>& models); - void registerCSTRModel(std::unordered_map>& models); + void registerInletModel(std::unordered_map>& models); + void registerOutletModel(std::unordered_map>& models); + + void registerGeneralRateModel(std::unordered_map>& models); + void registerLumpedRateModelWithPores(std::unordered_map>& models); + void registerLumpedRateModelWithoutPores(std::unordered_map>& models); + void registerCSTRModel(std::unordered_map>& models); #ifdef ENABLE_GRM_2D - void registerGeneralRateModel2D(std::unordered_map>& models); + void registerGeneralRateModel2D(std::unordered_map>& models); #endif namespace inlet @@ -63,12 +59,8 @@ namespace cadet model::registerInletModel(_modelCreators); model::registerOutletModel(_modelCreators); model::registerGeneralRateModel(_modelCreators); - model::registerGeneralRateModelDG(_modelCreators); - model::registerGeneralRateModelDGFV(_modelCreators); model::registerLumpedRateModelWithPores(_modelCreators); - model::registerLumpedRateModelWithPoresDG(_modelCreators); model::registerLumpedRateModelWithoutPores(_modelCreators); - model::registerLumpedRateModelWithoutPoresDG(_modelCreators); model::registerCSTRModel(_modelCreators); #ifdef ENABLE_GRM_2D @@ -92,7 +84,7 @@ namespace cadet template void ModelBuilder::registerModel(const std::string& name) { - _modelCreators[name] = [](UnitOpIdx uoId) { return new UnitOpModel_t(uoId); }; + _modelCreators[name] = [](UnitOpIdx uoId, IParameterProvider&) { return new UnitOpModel_t(uoId); }; } template @@ -188,7 +180,11 @@ namespace cadet // Call factory function (thanks to type erasure of std::function we can store // all factory functions in one container) - IUnitOperation* const model = it->second(uoId); + IUnitOperation* const model = it->second(uoId, paramProvider); + if (!model) { + LOG(Error) << "Failed to create unit type " << uoType << " for unit " << uoId; + return nullptr; + } if (!model->configureModelDiscretization(paramProvider, *this) || !model->configure(paramProvider)) { @@ -200,19 +196,19 @@ namespace cadet return model; } - IModel* ModelBuilder::createUnitOperation(const std::string& uoType, UnitOpIdx uoId) - { - const auto it = _modelCreators.find(uoType); - if (it == _modelCreators.end()) - { - // Model was not found - LOG(Error) << "Unknown unit type " << uoType << " for unit " << uoId; - return nullptr; - } - - IUnitOperation* const model = it->second(uoId); - return model; - } + //IModel* ModelBuilder::createUnitOperation(const std::string& uoType, UnitOpIdx uoId) + //{ + // const auto it = _modelCreators.find(uoType); + // if (it == _modelCreators.end()) + // { + // // Model was not found + // LOG(Error) << "Unknown unit type " << uoType << " for unit " << uoId; + // return nullptr; + // } + + // IUnitOperation* const model = it->second(uoId); + // return model; + //} void ModelBuilder::destroyUnitOperation(IModel* unitOp) { diff --git a/src/libcadet/ModelBuilderImpl.hpp b/src/libcadet/ModelBuilderImpl.hpp index e0f28a765..3b57706e6 100644 --- a/src/libcadet/ModelBuilderImpl.hpp +++ b/src/libcadet/ModelBuilderImpl.hpp @@ -52,7 +52,7 @@ class ModelBuilder : public IModelBuilder, public IConfigHelper virtual void destroySystem(IModelSystem* sys); virtual IModel* createUnitOperation(IParameterProvider& paramProvider, UnitOpIdx uoId); - virtual IModel* createUnitOperation(const std::string& uoType, UnitOpIdx uoId); + //virtual IModel* createUnitOperation(const std::string& uoType, UnitOpIdx uoId); virtual void destroyUnitOperation(IModel* unitOp); virtual void registerInletType(const std::string& name, std::function factory); @@ -91,7 +91,7 @@ class ModelBuilder : public IModelBuilder, public IConfigHelper ReactionModelFactory _reactionModels; //!< Factory for IDynamicReactionModel implementations ParameterDependenceFactory _paramDeps; //!< Factory for IParameterStateDependence implementations - typedef std::unordered_map> ModelFactoryContainer_t; + typedef std::unordered_map> ModelFactoryContainer_t; typedef std::unordered_map> InletFactoryContainer_t; typedef std::unordered_map> ExternalFunctionFactoryContainer_t; diff --git a/src/libcadet/model/GeneralRateModel.cpp b/src/libcadet/model/GeneralRateModel.cpp index 8e4881d45..595f8608f 100644 --- a/src/libcadet/model/GeneralRateModel.cpp +++ b/src/libcadet/model/GeneralRateModel.cpp @@ -3180,16 +3180,18 @@ namespace model template class GeneralRateModel; template class GeneralRateModel; -void registerGeneralRateModel(std::unordered_map>& models) +IUnitOperation* createAxialFVGRM(UnitOpIdx uoId) { typedef GeneralRateModel AxialGRM; - typedef GeneralRateModel RadialGRM; - models[AxialGRM::identifier()] = [](UnitOpIdx uoId) { return new AxialGRM(uoId); }; - models["GRM"] = [](UnitOpIdx uoId) { return new AxialGRM(uoId); }; + return new AxialGRM(uoId); +} + +IUnitOperation* createRadialFVGRM(UnitOpIdx uoId) +{ + typedef GeneralRateModel RadialGRM; - models[RadialGRM::identifier()] = [](UnitOpIdx uoId) { return new RadialGRM(uoId); }; - models["RGRM"] = [](UnitOpIdx uoId) { return new RadialGRM(uoId); }; + return new RadialGRM(uoId); } } // namespace model diff --git a/src/libcadet/model/GeneralRateModel.hpp b/src/libcadet/model/GeneralRateModel.hpp index ccb14bb95..c2c556adb 100644 --- a/src/libcadet/model/GeneralRateModel.hpp +++ b/src/libcadet/model/GeneralRateModel.hpp @@ -527,6 +527,9 @@ class GeneralRateModel : public UnitOperationBase extern template class GeneralRateModel; extern template class GeneralRateModel; +IUnitOperation* createAxialFVGRM(UnitOpIdx uoId); +IUnitOperation* createRadialFVGRM(UnitOpIdx uoId); + } // namespace model } // namespace cadet diff --git a/src/libcadet/model/GeneralRateModel2D.cpp b/src/libcadet/model/GeneralRateModel2D.cpp index 01d201770..7304f4235 100644 --- a/src/libcadet/model/GeneralRateModel2D.cpp +++ b/src/libcadet/model/GeneralRateModel2D.cpp @@ -2568,10 +2568,10 @@ int GeneralRateModel2D::Exporter::writeOutlet(double* buffer) const } -void registerGeneralRateModel2D(std::unordered_map>& models) +void registerGeneralRateModel2D(std::unordered_map>& models) { - models[GeneralRateModel2D::identifier()] = [](UnitOpIdx uoId) { return new GeneralRateModel2D(uoId); }; - models["GRM2D"] = [](UnitOpIdx uoId) { return new GeneralRateModel2D(uoId); }; + models[GeneralRateModel2D::identifier()] = [](UnitOpIdx uoId, IParameterProvider&) { return new GeneralRateModel2D(uoId); }; + models["GRM2D"] = [](UnitOpIdx uoId, IParameterProvider&) { return new GeneralRateModel2D(uoId); }; } } // namespace model diff --git a/src/libcadet/model/GeneralRateModelBuilder.cpp b/src/libcadet/model/GeneralRateModelBuilder.cpp new file mode 100644 index 000000000..4d091d568 --- /dev/null +++ b/src/libcadet/model/GeneralRateModelBuilder.cpp @@ -0,0 +1,91 @@ +#include "model/GeneralRateModel.hpp" +#include "model/GeneralRateModelDG.hpp" +#include "model/GeneralRateModelDGFV.hpp" +#include "LoggingUtils.hpp" +#include "Logging.hpp" + + +namespace cadet +{ + namespace model + { + + IUnitOperation* selectAxialFlowDiscretizationGRM(UnitOpIdx uoId, IParameterProvider& paramProvider) + { + IUnitOperation* model = nullptr; + + paramProvider.pushScope("discretization"); + + if (paramProvider.exists("SPATIAL_METHOD")) { + + const std::string discName = paramProvider.getString("SPATIAL_METHOD"); + + if (discName == "DG") + model = new GeneralRateModelDG(uoId); + else if (discName == "FV") + model = createAxialFVGRM(uoId); + else + { + LOG(Error) << "Unknown discretization type " << discName << " for unit " << uoId; + } + } + else { + model = createAxialFVGRM(uoId); + } + + paramProvider.popScope(); + + return model; + } + + IUnitOperation* selectRadialFlowDiscretizationGRM(UnitOpIdx uoId, IParameterProvider& paramProvider) + { + IUnitOperation* model = nullptr; + + paramProvider.pushScope("discretization"); + + if (paramProvider.exists("SPATIAL_METHOD")) { + + const std::string discName = paramProvider.getString("SPATIAL_METHOD"); + + if (discName == "DG") + { + LOG(Error) << "Radial flow not implemented for DG discretization yet, was called for unit " << uoId; + } + else if (discName == "FV") + model = createRadialFVGRM(uoId); + else + { + LOG(Error) << "Unknown discretization type " << discName << " for unit " << uoId; + } + } + else { + model = createRadialFVGRM(uoId); + } + + paramProvider.popScope(); + + return model; + } + + void registerGeneralRateModel(std::unordered_map>& models) + { + typedef GeneralRateModel AxialGRM; + typedef GeneralRateModel RadialGRM; + + models[GeneralRateModelDG::identifier()] = selectAxialFlowDiscretizationGRM; + models["GRM_DG"] = selectAxialFlowDiscretizationGRM; + + models[GeneralRateModelDGFV::identifier()] = selectAxialFlowDiscretizationGRM; + models["GRM_DGFV"] = selectAxialFlowDiscretizationGRM; + + models[AxialGRM::identifier()] = selectAxialFlowDiscretizationGRM; + models["GRM"] = selectAxialFlowDiscretizationGRM; + + models[RadialGRM::identifier()] = selectRadialFlowDiscretizationGRM; + models["RGRM"] = selectRadialFlowDiscretizationGRM; + } + + } // namespace model + +} // namespace cadet \ No newline at end of file diff --git a/src/libcadet/model/GeneralRateModelDG.cpp b/src/libcadet/model/GeneralRateModelDG.cpp index 75430a360..d6b0d88fc 100644 --- a/src/libcadet/model/GeneralRateModelDG.cpp +++ b/src/libcadet/model/GeneralRateModelDG.cpp @@ -134,95 +134,103 @@ bool GeneralRateModelDG::configureModelDiscretization(IParameterProvider& paramP paramProvider.pushScope("discretization"); - _disc.nCol = paramProvider.getInt("NCOL"); - if (_disc.nCol < 1) - throw InvalidParameterException("Number of column cells must be at least 1!"); - + if (paramProvider.exists("POLYDEG")) + _disc.polyDeg = paramProvider.getInt("POLYDEG"); + else + _disc.polyDeg = 4u; // default value if (paramProvider.getInt("POLYDEG") < 1) throw InvalidParameterException("Polynomial degree must be at least 1!"); + else if (_disc.polyDeg < 3) + LOG(Warning) << "Polynomial degree > 2 in bulk discretization (cf. POLYDEG) is always recommended for performance reasons."; + + _disc.nNodes = _disc.polyDeg + 1; + + if (paramProvider.exists("NELEM")) + _disc.nCol = paramProvider.getInt("NELEM"); + else if (paramProvider.exists("NCOL")) + _disc.nCol = std::max(1u, paramProvider.getInt("NCOL") / _disc.nNodes); // number of elements is rounded down else - _disc.polyDeg = paramProvider.getInt("POLYDEG"); + throw InvalidParameterException("Specify field NELEM (or NCOL)"); + + if (_disc.nCol < 1) + throw InvalidParameterException("Number of column elements must be at least 1!"); - if (_disc.polyDeg < 3) - LOG(Warning) << "Polynomial degree > 2 in bulk discretization (cf. POLYDEG) is always recommended for performance reasons."; - _disc.nNodes = _disc.polyDeg + 1u; _disc.nPoints = _disc.nNodes * _disc.nCol; - const int polynomial_integration_mode = paramProvider.getInt("EXACT_INTEGRATION"); - _disc.exactInt = static_cast(polynomial_integration_mode); // only integration mode 0 applies the inexact collocated diagonal LGL mass matrix - const std::vector nParCell = paramProvider.getIntArray("NPARCELL"); - const std::vector parPolyDeg = paramProvider.getIntArray("PARPOLYDEG"); - const std::vector parExactInt = paramProvider.getBoolArray("PAR_EXACT_INTEGRATION"); + int polynomial_integration_mode = 0; + if (paramProvider.exists("EXACT_INTEGRATION")) + polynomial_integration_mode = paramProvider.getInt("EXACT_INTEGRATION"); + _disc.exactInt = static_cast(polynomial_integration_mode); // only integration mode 0 applies the inexact collocated diagonal LGL mass matrix const std::vector nBound = paramProvider.getIntArray("NBOUND"); if (nBound.size() < _disc.nComp) throw InvalidParameterException("Field NBOUND contains too few elements (NCOMP = " + std::to_string(_disc.nComp) + " required)"); + if ((nBound.size() > _disc.nComp) && (nBound.size() < _disc.nComp * _disc.nParType)) + throw InvalidParameterException("Field NBOUND must have NCOMP (" + std::to_string(_disc.nComp) + ") or NCOMP * NPARTYPE (" + std::to_string(_disc.nComp * _disc.nParType) + ") entries"); if (paramProvider.exists("NPARTYPE")) _disc.nParType = paramProvider.getInt("NPARTYPE"); - else - { - // Infer number of particle types - _disc.nParType = std::max({ nBound.size() / _disc.nComp, nParCell.size(), parPolyDeg.size(), parExactInt.size() }); - } - - if ((parExactInt.size() > 1) && (parExactInt.size() < _disc.nParType)) - throw InvalidParameterException("Field PAR_EXACT_INTEGRATION must have 1 or NPARTYPE (" + std::to_string(_disc.nParType) + ") entries"); - if ((nParCell.size() > 1) && (nParCell.size() < _disc.nParType)) - throw InvalidParameterException("Field NPARCELL must have 1 or NPARTYPE (" + std::to_string(_disc.nParType) + ") entries"); - if ((parPolyDeg.size() > 1) && (parPolyDeg.size() < _disc.nParType)) - throw InvalidParameterException("Field PARPOLYDEG must have 1 or NPARTYPE (" + std::to_string(_disc.nParType) + ") entries"); - - _disc.nParCell = new unsigned int[_disc.nParType]; - if (nParCell.size() < _disc.nParType) - { - // Multiplex number of particle cells to all particle types - for (unsigned int i = 0; i < _disc.nParType; ++i) - std::fill(_disc.nParCell, _disc.nParCell + _disc.nParType, nParCell[0]); - } - else - std::copy_n(nParCell.begin(), _disc.nParType, _disc.nParCell); - - _disc.parPolyDeg = new unsigned int[_disc.nParType]; - if (parPolyDeg.size() < _disc.nParType) - { - if (parPolyDeg[0] < 1) { - throw InvalidParameterException("Particle polynomial degree must be at least 1!"); - } - else { - // Multiplex polynomial degree of particle elements to all particle types + else // Infer number of particle types + _disc.nParType = nBound.size() / _disc.nComp; + + std::vector parPolyDeg(_disc.nParType); + std::vector ParNelem(_disc.nParType); + std::vector parExactInt(_disc.nParType); + + if (paramProvider.exists("PAR_POLYDEG")) + { + parPolyDeg = paramProvider.getIntArray("PAR_POLYDEG"); + if ((std::any_of(parPolyDeg.begin(), parPolyDeg.end(), [](int value) { return value < 1; }))) + throw InvalidParameterException("Particle polynomial degrees must be at least 1!"); + ParNelem = paramProvider.getIntArray("PAR_NELEM"); + if ((std::any_of(ParNelem.begin(), ParNelem.end(), [](int value) { return value < 1; }))) + throw InvalidParameterException("Particle number of elements must be at least 1!"); + parExactInt = paramProvider.getBoolArray("PAR_EXACT_INTEGRATION"); + if ((std::any_of(parExactInt.begin(), parExactInt.end(), [](bool value) { return !value; }))) + LOG(Warning) << "Inexact integration method (cf. PAR_EXACT_INTEGRATION) in particles might add severe! stiffness to the system and disables consistent initialization!"; + + if (parPolyDeg.size() == 1) + { + // Multiplex number of particle shells to all particle types for (unsigned int i = 0; i < _disc.nParType; ++i) std::fill(_disc.parPolyDeg, _disc.parPolyDeg + _disc.nParType, parPolyDeg[0]); - } - } - else { - for (unsigned int parType = 0; parType < _disc.nParType; parType++) + else if (parPolyDeg.size() < _disc.nParType) + throw InvalidParameterException("Field PARPOLYDEG must have 1 or NPARTYPE (" + std::to_string(_disc.nParType) + ") entries"); + else + std::copy_n(parPolyDeg.begin(), _disc.nParType, _disc.parPolyDeg); + if (ParNelem.size() == 1) { - if (parPolyDeg[parType] < 1) - throw InvalidParameterException("Particle polynomial degree(s) must be at least 1!"); - if (_disc.nParCell[parType] < 1) - throw InvalidParameterException("Number of particle cell(s) must be at least 1!"); - - _disc.parPolyDeg[parType] = parPolyDeg[parType]; + // Multiplex number of particle shells to all particle types + for (unsigned int i = 0; i < _disc.nParType; ++i) + std::fill(_disc.nParCell, _disc.nParCell + _disc.nParType, ParNelem[0]); } + else if (ParNelem.size() < _disc.nParType) + throw InvalidParameterException("Field PAR_NELEM must have 1 or NPARTYPE (" + std::to_string(_disc.nParType) + ") entries"); + else + std::copy_n(ParNelem.begin(), _disc.nParType, _disc.nParCell); + if (parExactInt.size() == 1) + { + // Multiplex number of particle shells to all particle types + for (unsigned int i = 0; i < _disc.nParType; ++i) + std::fill(_disc.parExactInt, _disc.parExactInt + _disc.nParType, parExactInt[0]); + } + else if (parExactInt.size() < _disc.nParType) + throw InvalidParameterException("Field PAR_EXACT_INTEGRATION must have 1 or NPARTYPE (" + std::to_string(_disc.nParType) + ") entries"); + else + std::copy_n(parExactInt.begin(), _disc.nParType, _disc.parExactInt); } - - _disc.parExactInt = new bool[_disc.nParType]; - if (parExactInt.size() < _disc.nParType) + else if (paramProvider.exists("NPAR")) { - // Multiplex exact/inexact integration of particle elements to all particle types - for (unsigned int i = 0; i < _disc.nParType; ++i) { - std::fill(_disc.parExactInt, _disc.parExactInt + _disc.nParType, parExactInt[0]); - if (!_disc.parExactInt[i]) - LOG(Warning) << "Inexact integration method (cf. PAR_EXACT_INTEGRATION) in particles might add severe! stiffness to the system and disables consistent initialization!"; - } + const std::vector nParPoints = paramProvider.getIntArray("NPAR"); + if ((nParPoints.size() > 1) && (nParPoints.size() < _disc.nParType)) + throw InvalidParameterException("Field NPAR must have 1 or NPARTYPE (" + std::to_string(_disc.nParType) + ") entries"); + + for (unsigned int par = 0; par < _disc.nParType; par++) + parPolyDeg[par] = std::max(1, std::min(nParPoints[par] - 1, 4)); } else - std::copy_n(parExactInt.begin(), _disc.nParType, _disc.parExactInt); - - if ((nBound.size() > _disc.nComp) && (nBound.size() < _disc.nComp * _disc.nParType)) - throw InvalidParameterException("Field NBOUND must have NCOMP (" + std::to_string(_disc.nComp) + ") or NCOMP * NPARTYPE (" + std::to_string(_disc.nComp * _disc.nParType) + ") entries"); + throw InvalidParameterException("Specify field PARPOLYDEG (or NPAR)"); // Compute discretization operators and initialize containers _disc.initializeDG(); @@ -2408,18 +2416,3 @@ int GeneralRateModelDG::Exporter::writeOutlet(double* buffer) const #include "model/GeneralRateModelDG-InitialConditions.cpp" #include "model/GeneralRateModelDG-LinearSolver.cpp" -namespace cadet -{ - -namespace model -{ - -void registerGeneralRateModelDG(std::unordered_map>& models) -{ - models[GeneralRateModelDG::identifier()] = [](UnitOpIdx uoId) { return new GeneralRateModelDG(uoId); }; - models["GRM_DG"] = [](UnitOpIdx uoId) { return new GeneralRateModelDG(uoId); }; -} - -} // namespace model - -} // namespace cadet diff --git a/src/libcadet/model/GeneralRateModelDG.hpp b/src/libcadet/model/GeneralRateModelDG.hpp index 35971edf5..027866d34 100644 --- a/src/libcadet/model/GeneralRateModelDG.hpp +++ b/src/libcadet/model/GeneralRateModelDG.hpp @@ -821,8 +821,8 @@ class GeneralRateModelDG : public UnitOperationBase virtual bool hasParticleMobilePhase() const CADET_NOEXCEPT { return true; } virtual bool hasSolidPhase() const CADET_NOEXCEPT { return _disc.strideBound[_disc.nParType] > 0; } virtual bool hasVolume() const CADET_NOEXCEPT { return false; } - virtual bool hasSmoothnessIndicator() const CADET_NOEXCEPT { return _model._convDispOp.hasSmoothnessIndicator(); } virtual bool isParticleLumped() const CADET_NOEXCEPT { return false; } + virtual bool hasSmoothnessIndicator() const CADET_NOEXCEPT { return _model._convDispOp.hasSmoothnessIndicator(); } virtual unsigned int primaryPolynomialDegree() const CADET_NOEXCEPT { return _disc.polyDeg; } virtual unsigned int secondaryPolynomialDegree() const CADET_NOEXCEPT { return 0; } diff --git a/src/libcadet/model/GeneralRateModelDGFV.cpp b/src/libcadet/model/GeneralRateModelDGFV.cpp index 799337400..0ae754c6d 100644 --- a/src/libcadet/model/GeneralRateModelDGFV.cpp +++ b/src/libcadet/model/GeneralRateModelDGFV.cpp @@ -3192,19 +3192,3 @@ int GeneralRateModelDGFV::Exporter::writeOutlet(double* buffer) const #include "model/GeneralRateModelDGFV-InitialConditions.cpp" #include "model/GeneralRateModelDGFV-LinearSolver.cpp" - -namespace cadet -{ - -namespace model -{ - -void registerGeneralRateModelDGFV(std::unordered_map>& models) -{ - models[GeneralRateModelDGFV::identifier()] = [](UnitOpIdx uoId) { return new GeneralRateModelDGFV(uoId); }; - models["GRM_DGFV"] = [](UnitOpIdx uoId) { return new GeneralRateModelDGFV(uoId); }; -} - -} // namespace model - -} // namespace cadet diff --git a/src/libcadet/model/GeneralRateModelDGFV.hpp b/src/libcadet/model/GeneralRateModelDGFV.hpp index dc1963069..45df14c5d 100644 --- a/src/libcadet/model/GeneralRateModelDGFV.hpp +++ b/src/libcadet/model/GeneralRateModelDGFV.hpp @@ -696,8 +696,8 @@ class GeneralRateModelDGFV : public UnitOperationBase virtual bool hasParticleMobilePhase() const CADET_NOEXCEPT { return true; } virtual bool hasSolidPhase() const CADET_NOEXCEPT { return _disc.strideBound[_disc.nParType] > 0; } virtual bool hasVolume() const CADET_NOEXCEPT { return false; } - virtual bool hasSmoothnessIndicator() const CADET_NOEXCEPT { return false; } virtual bool isParticleLumped() const CADET_NOEXCEPT { return false; } + virtual bool hasSmoothnessIndicator() const CADET_NOEXCEPT { return false; } virtual unsigned int primaryPolynomialDegree() const CADET_NOEXCEPT { return _disc.polyDeg; } virtual unsigned int secondaryPolynomialDegree() const CADET_NOEXCEPT { return 0; } diff --git a/src/libcadet/model/InletModel.cpp b/src/libcadet/model/InletModel.cpp index 3bb9908f0..588aeeb34 100644 --- a/src/libcadet/model/InletModel.cpp +++ b/src/libcadet/model/InletModel.cpp @@ -444,9 +444,9 @@ int InletModel::Exporter::writeOutlet(double* buffer) const } -void registerInletModel(std::unordered_map>& models) +void registerInletModel(std::unordered_map>& models) { - models[InletModel::identifier()] = [](UnitOpIdx uoId) { return new InletModel(uoId); }; + models[InletModel::identifier()] = [](UnitOpIdx uoId, IParameterProvider&) { return new InletModel(uoId); }; } } // namespace model diff --git a/src/libcadet/model/LumpedRateModelWithPores.cpp b/src/libcadet/model/LumpedRateModelWithPores.cpp index 463adaf01..47d7f9cb2 100644 --- a/src/libcadet/model/LumpedRateModelWithPores.cpp +++ b/src/libcadet/model/LumpedRateModelWithPores.cpp @@ -1762,16 +1762,18 @@ namespace model template class LumpedRateModelWithPores; template class LumpedRateModelWithPores; -void registerLumpedRateModelWithPores(std::unordered_map>& models) +IUnitOperation* createAxialFVLRMP(UnitOpIdx uoId) { typedef LumpedRateModelWithPores AxialLRMP; - typedef LumpedRateModelWithPores RadialLRMP; - models[AxialLRMP::identifier()] = [](UnitOpIdx uoId) { return new AxialLRMP(uoId); }; - models["LRMP"] = [](UnitOpIdx uoId) { return new AxialLRMP(uoId); }; + return new AxialLRMP(uoId); +} + +IUnitOperation* createRadialFVLRMP(UnitOpIdx uoId) +{ + typedef LumpedRateModelWithPores RadialLRMP; - models[RadialLRMP::identifier()] = [](UnitOpIdx uoId) { return new RadialLRMP(uoId); }; - models["RLRMP"] = [](UnitOpIdx uoId) { return new RadialLRMP(uoId); }; + return new RadialLRMP(uoId); } } // namespace model diff --git a/src/libcadet/model/LumpedRateModelWithPores.hpp b/src/libcadet/model/LumpedRateModelWithPores.hpp index 3f09858f0..839171748 100644 --- a/src/libcadet/model/LumpedRateModelWithPores.hpp +++ b/src/libcadet/model/LumpedRateModelWithPores.hpp @@ -448,6 +448,9 @@ class LumpedRateModelWithPores : public UnitOperationBase extern template class LumpedRateModelWithPores; extern template class LumpedRateModelWithPores; +IUnitOperation* createAxialFVLRMP(UnitOpIdx uoId); +IUnitOperation* createRadialFVLRMP(UnitOpIdx uoId); + } // namespace model } // namespace cadet diff --git a/src/libcadet/model/LumpedRateModelWithPoresBuilder.cpp b/src/libcadet/model/LumpedRateModelWithPoresBuilder.cpp new file mode 100644 index 000000000..2b9c459c4 --- /dev/null +++ b/src/libcadet/model/LumpedRateModelWithPoresBuilder.cpp @@ -0,0 +1,87 @@ +#include "model/LumpedRateModelWithPores.hpp" +#include "model/LumpedRateModelWithPoresDG.hpp" +#include "LoggingUtils.hpp" +#include "Logging.hpp" + + +namespace cadet +{ + namespace model + { + + IUnitOperation* selectAxialFlowDiscretizationLRMP(UnitOpIdx uoId, IParameterProvider& paramProvider) + { + IUnitOperation* model = nullptr; + + paramProvider.pushScope("discretization"); + + if (paramProvider.exists("SPATIAL_METHOD")) { + + const std::string discName = paramProvider.getString("SPATIAL_METHOD"); + + if (discName == "DG") + model = new LumpedRateModelWithPoresDG(uoId); + else if (discName == "FV") + model = createAxialFVLRMP(uoId); + else + { + LOG(Error) << "Unknown discretization type " << discName << " for unit " << uoId; + } + } + else { + model = createAxialFVLRMP(uoId); + } + + paramProvider.popScope(); + + return model; + } + + IUnitOperation* selectRadialFlowDiscretizationLRMP(UnitOpIdx uoId, IParameterProvider& paramProvider) + { + IUnitOperation* model = nullptr; + + paramProvider.pushScope("discretization"); + + if (paramProvider.exists("SPATIAL_METHOD")) { + + const std::string discName = paramProvider.getString("SPATIAL_METHOD"); + + if (discName == "DG") + { + LOG(Error) << "Radial flow not implemented for DG discretization yet, was called for unit " << uoId; + } + else if (discName == "FV") + model = createRadialFVLRMP(uoId); + else + { + LOG(Error) << "Unknown discretization type " << discName << " for unit " << uoId; + } + } + else { + model = createRadialFVLRMP(uoId); + } + + paramProvider.popScope(); + + return model; + } + + void registerLumpedRateModelWithPores(std::unordered_map>& models) + { + typedef LumpedRateModelWithPores AxialLRMP; + typedef LumpedRateModelWithPores RadialLRMP; + + models[LumpedRateModelWithPoresDG::identifier()] = selectAxialFlowDiscretizationLRMP; + models["LRMP_DG"] = selectAxialFlowDiscretizationLRMP; + + models[AxialLRMP::identifier()] = selectAxialFlowDiscretizationLRMP; + models["LRMP"] = selectAxialFlowDiscretizationLRMP; + + models[RadialLRMP::identifier()] = selectRadialFlowDiscretizationLRMP; + models["RLRMP"] = selectRadialFlowDiscretizationLRMP; + } + + } // namespace model + +} // namespace cadet \ No newline at end of file diff --git a/src/libcadet/model/LumpedRateModelWithPoresDG.cpp b/src/libcadet/model/LumpedRateModelWithPoresDG.cpp index 1d722bd0b..482856b79 100644 --- a/src/libcadet/model/LumpedRateModelWithPoresDG.cpp +++ b/src/libcadet/model/LumpedRateModelWithPoresDG.cpp @@ -103,18 +103,32 @@ bool LumpedRateModelWithPoresDG::configureModelDiscretization(IParameterProvider paramProvider.pushScope("discretization"); - if (paramProvider.getInt("NCOL") < 1) - throw InvalidParameterException("Number of column cells must be at least 1!"); - _disc.nCol = paramProvider.getInt("NCOL"); - + if (paramProvider.exists("POLYDEG")) + _disc.polyDeg = paramProvider.getInt("POLYDEG"); + else + _disc.polyDeg = 4u; // default value if (paramProvider.getInt("POLYDEG") < 1) throw InvalidParameterException("Polynomial degree must be at least 1!"); - else - _disc.polyDeg = paramProvider.getInt("POLYDEG"); + else if (_disc.polyDeg < 3) + LOG(Warning) << "Polynomial degree > 2 in bulk discretization (cf. POLYDEG) is always recommended for performance reasons."; + _disc.nNodes = _disc.polyDeg + 1; + + if (paramProvider.exists("NELEM")) + _disc.nCol = paramProvider.getInt("NELEM"); + else if (paramProvider.exists("NCOL")) + _disc.nCol = std::max(1u, paramProvider.getInt("NCOL") / _disc.nNodes); // number of elements is rounded down + else + throw InvalidParameterException("Specify field NELEM (or NCOL)"); + + if (_disc.nCol < 1) + throw InvalidParameterException("Number of column elements must be at least 1!"); + _disc.nPoints = _disc.nNodes * _disc.nCol; - const int polynomial_integration_mode = paramProvider.getInt("EXACT_INTEGRATION"); + int polynomial_integration_mode = 0; + if(paramProvider.exists("EXACT_INTEGRATION")) + polynomial_integration_mode = paramProvider.getInt("EXACT_INTEGRATION"); _disc.exactInt = static_cast(polynomial_integration_mode); // only integration mode 0 applies the inexact collocated diagonal LGL mass matrix const std::vector nBound = paramProvider.getIntArray("NBOUND"); @@ -1509,15 +1523,15 @@ int LumpedRateModelWithPoresDG::Exporter::writeOutlet(double* buffer) const namespace cadet { -namespace model -{ + namespace model + { -void registerLumpedRateModelWithPoresDG(std::unordered_map>& models) -{ - models[LumpedRateModelWithPoresDG::identifier()] = [](UnitOpIdx uoId) { return new LumpedRateModelWithPoresDG(uoId); }; - models["LRMPDG"] = [](UnitOpIdx uoId) { return new LumpedRateModelWithPoresDG(uoId); }; -} + void registerLumpedRateModelWithPoresDG(std::unordered_map>& models) + { + models[LumpedRateModelWithPoresDG::identifier()] = [](UnitOpIdx uoId) { return new LumpedRateModelWithPoresDG(uoId); }; + models["LRMPDG"] = [](UnitOpIdx uoId) { return new LumpedRateModelWithPoresDG(uoId); }; + } -} // namespace model + } // namespace model } // namespace cadet diff --git a/src/libcadet/model/LumpedRateModelWithoutPores.cpp b/src/libcadet/model/LumpedRateModelWithoutPores.cpp index 0e981b78e..fd81baeae 100644 --- a/src/libcadet/model/LumpedRateModelWithoutPores.cpp +++ b/src/libcadet/model/LumpedRateModelWithoutPores.cpp @@ -1911,20 +1911,20 @@ int LumpedRateModelWithoutPores::Exporter::writeOutlet(double* return _disc.nComp; } - -void registerLumpedRateModelWithoutPores(std::unordered_map>& models) +IUnitOperation* createAxialFVLRM(UnitOpIdx uoId) { typedef LumpedRateModelWithoutPores AxialLRM; - typedef LumpedRateModelWithoutPores RadialLRM; - models[AxialLRM::identifier()] = [](UnitOpIdx uoId) { return new AxialLRM(uoId); }; - models["LRM"] = [](UnitOpIdx uoId) { return new AxialLRM(uoId); }; - models["DPFR"] = [](UnitOpIdx uoId) { return new AxialLRM(uoId); }; + return new AxialLRM(uoId); +} + +IUnitOperation* createRadialFVLRM(UnitOpIdx uoId) +{ + typedef LumpedRateModelWithoutPores RadialLRM; - models[RadialLRM::identifier()] = [](UnitOpIdx uoId) { return new RadialLRM(uoId); }; - models["RLRM"] = [](UnitOpIdx uoId) { return new RadialLRM(uoId); }; + return new RadialLRM(uoId); } } // namespace model -} // namespace cadet +} // namespace cadet \ No newline at end of file diff --git a/src/libcadet/model/LumpedRateModelWithoutPores.hpp b/src/libcadet/model/LumpedRateModelWithoutPores.hpp index b58d32efd..3567dcf0c 100644 --- a/src/libcadet/model/LumpedRateModelWithoutPores.hpp +++ b/src/libcadet/model/LumpedRateModelWithoutPores.hpp @@ -347,6 +347,9 @@ class LumpedRateModelWithoutPores : public UnitOperationBase }; }; +IUnitOperation* createAxialFVLRM(UnitOpIdx uoId); +IUnitOperation* createRadialFVLRM(UnitOpIdx uoId); + } // namespace model } // namespace cadet diff --git a/src/libcadet/model/LumpedRateModelWithoutPoresBuilder.cpp b/src/libcadet/model/LumpedRateModelWithoutPoresBuilder.cpp new file mode 100644 index 000000000..0b5a9ad8c --- /dev/null +++ b/src/libcadet/model/LumpedRateModelWithoutPoresBuilder.cpp @@ -0,0 +1,88 @@ +#include "model/LumpedRateModelWithoutPores.hpp" +#include "model/LumpedRateModelWithoutPoresDG.hpp" +#include "LoggingUtils.hpp" +#include "Logging.hpp" + + +namespace cadet +{ +namespace model +{ + + IUnitOperation* selectAxialFlowDiscretizationLRM(UnitOpIdx uoId, IParameterProvider& paramProvider) + { + IUnitOperation* model = nullptr; + + paramProvider.pushScope("discretization"); + + if (paramProvider.exists("SPATIAL_METHOD")) { + + const std::string discName = paramProvider.getString("SPATIAL_METHOD"); + + if(discName == "DG") + model = new LumpedRateModelWithoutPoresDG(uoId); + else if (discName == "FV") + model = createAxialFVLRM(uoId); + else + { + LOG(Error) << "Unknown discretization type " << discName << " for unit " << uoId; + } + } + else { + model = createAxialFVLRM(uoId); + } + + paramProvider.popScope(); + + return model; + } + + IUnitOperation* selectRadialFlowDiscretizationLRM(UnitOpIdx uoId, IParameterProvider& paramProvider) + { + IUnitOperation* model = nullptr; + + paramProvider.pushScope("discretization"); + + if (paramProvider.exists("SPATIAL_METHOD")) { + + const std::string discName = paramProvider.getString("SPATIAL_METHOD"); + + if (discName == "DG") + { + LOG(Error) << "Radial flow not implemented for DG discretization yet, was called for unit " << uoId; + } + else if (discName == "FV") + model = createRadialFVLRM(uoId); + else + { + LOG(Error) << "Unknown discretization type " << discName << " for unit " << uoId; + } + } + else { + model = createRadialFVLRM(uoId); + } + + paramProvider.popScope(); + + return model; + } + +void registerLumpedRateModelWithoutPores(std::unordered_map>& models) +{ + typedef LumpedRateModelWithoutPores AxialLRM; + typedef LumpedRateModelWithoutPores RadialLRM; + + models[LumpedRateModelWithoutPoresDG::identifier()] = selectAxialFlowDiscretizationLRM; + models["LRM_DG"] = selectAxialFlowDiscretizationLRM; + + models[AxialLRM::identifier()] = selectAxialFlowDiscretizationLRM; + models["LRM"] = selectAxialFlowDiscretizationLRM; + models["DPFR"] = selectAxialFlowDiscretizationLRM; + + models[RadialLRM::identifier()] = selectRadialFlowDiscretizationLRM; + models["RLRM"] = selectRadialFlowDiscretizationLRM; +} + +} // namespace model + +} // namespace cadet \ No newline at end of file diff --git a/src/libcadet/model/LumpedRateModelWithoutPoresDG.cpp b/src/libcadet/model/LumpedRateModelWithoutPoresDG.cpp index 258531189..5c5090a66 100644 --- a/src/libcadet/model/LumpedRateModelWithoutPoresDG.cpp +++ b/src/libcadet/model/LumpedRateModelWithoutPoresDG.cpp @@ -102,20 +102,34 @@ namespace cadet paramProvider.pushScope("discretization"); - const int polynomial_integration_mode = paramProvider.getInt("EXACT_INTEGRATION"); - _disc.exactInt = static_cast(polynomial_integration_mode); // only integration mode 0 applies the inexact collocated diagonal LGL mass matrix - + if (paramProvider.exists("POLYDEG")) + _disc.polyDeg = paramProvider.getInt("POLYDEG"); + else + _disc.polyDeg = 4u; // default value if (paramProvider.getInt("POLYDEG") < 1) throw InvalidParameterException("Polynomial degree must be at least 1!"); - _disc.polyDeg = paramProvider.getInt("POLYDEG"); - _disc.nNodes = _disc.polyDeg + 1u; + else if (_disc.polyDeg < 3) + LOG(Warning) << "Polynomial degree > 2 in bulk discretization (cf. POLYDEG) is always recommended for performance reasons."; + + _disc.nNodes = _disc.polyDeg + 1; + + if (paramProvider.exists("NELEM")) + _disc.nCol = paramProvider.getInt("NELEM"); + else if (paramProvider.exists("NCOL")) + _disc.nCol = std::max(1u, paramProvider.getInt("NCOL") / _disc.nNodes); // number of elements is rounded down + else + throw InvalidParameterException("Specify field NELEM (or NCOL)"); - if (paramProvider.getInt("NCOL") < 1) - throw InvalidParameterException("Number of column cells must be at least 1!"); - _disc.nCol = paramProvider.getInt("NCOL"); + if (_disc.nCol < 1) + throw InvalidParameterException("Number of column elements must be at least 1!"); _disc.nPoints = _disc.nNodes * _disc.nCol; + int polynomial_integration_mode = 0; + if (paramProvider.exists("EXACT_INTEGRATION")) + polynomial_integration_mode = paramProvider.getInt("EXACT_INTEGRATION"); + _disc.exactInt = static_cast(polynomial_integration_mode); // only integration mode 0 applies the inexact collocated diagonal LGL mass matrix + const std::vector nBound = paramProvider.getIntArray("NBOUND"); if (nBound.size() < _disc.nComp) throw InvalidParameterException("Field NBOUND contains too few elements (NCOMP = " + std::to_string(_disc.nComp) + " required)"); @@ -1888,13 +1902,6 @@ namespace cadet } - void registerLumpedRateModelWithoutPoresDG(std::unordered_map>& models) - { - models[LumpedRateModelWithoutPoresDG::identifier()] = [](UnitOpIdx uoId) { return new LumpedRateModelWithoutPoresDG(uoId); }; - models["LRM_DG"] = [](UnitOpIdx uoId) { return new LumpedRateModelWithoutPoresDG(uoId); }; - models["DPFR_DG"] = [](UnitOpIdx uoId) { return new LumpedRateModelWithoutPoresDG(uoId); }; - } - } // namespace model } // namespace cadet diff --git a/src/libcadet/model/OutletModel.cpp b/src/libcadet/model/OutletModel.cpp index 7354148cf..1f6c4fcf2 100644 --- a/src/libcadet/model/OutletModel.cpp +++ b/src/libcadet/model/OutletModel.cpp @@ -256,9 +256,9 @@ int OutletModel::Exporter::writeOutlet(double* buffer) const } -void registerOutletModel(std::unordered_map>& models) +void registerOutletModel(std::unordered_map>& models) { - models[OutletModel::identifier()] = [](UnitOpIdx uoId) { return new OutletModel(uoId); }; + models[OutletModel::identifier()] = [](UnitOpIdx uoId, IParameterProvider&) { return new OutletModel(uoId); }; } } // namespace model diff --git a/src/libcadet/model/StirredTankModel.cpp b/src/libcadet/model/StirredTankModel.cpp index b73ae97b8..059d7e049 100644 --- a/src/libcadet/model/StirredTankModel.cpp +++ b/src/libcadet/model/StirredTankModel.cpp @@ -2023,9 +2023,9 @@ int CSTRModel::Exporter::writeOutlet(double* buffer) const -void registerCSTRModel(std::unordered_map>& models) +void registerCSTRModel(std::unordered_map>& models) { - models[CSTRModel::identifier()] = [](UnitOpIdx uoId) { return new CSTRModel(uoId); }; + models[CSTRModel::identifier()] = [](UnitOpIdx uoId, IParameterProvider&) { return new CSTRModel(uoId); }; } } // namespace model diff --git a/test/CSTR-Residual.cpp b/test/CSTR-Residual.cpp index 1421a774d..55bd296ea 100644 --- a/test/CSTR-Residual.cpp +++ b/test/CSTR-Residual.cpp @@ -40,7 +40,7 @@ cadet::model::CSTRModel* createAndConfigureCSTR(cadet::IModelBuilder& mb, cadet::JsonParameterProvider& jpp) { // Create a CSTR - cadet::IModel* const iCstr = mb.createUnitOperation("CSTR", 0); + cadet::IModel* const iCstr = mb.createUnitOperation(jpp, 0); REQUIRE(nullptr != iCstr); cadet::model::CSTRModel* const cstr = reinterpret_cast(iCstr); diff --git a/test/ColumnTests.cpp b/test/ColumnTests.cpp index 38d761fc1..796b79224 100644 --- a/test/ColumnTests.cpp +++ b/test/ColumnTests.cpp @@ -56,7 +56,7 @@ namespace inline cadet::IUnitOperation* createAndConfigureUnit(const std::string& uoType, cadet::IModelBuilder& mb, cadet::JsonParameterProvider& jpp, int wenoOrder) { // Create a unit - cadet::IModel* const iUnit = mb.createUnitOperation(uoType, 0); + cadet::IModel* const iUnit = mb.createUnitOperation(jpp, 0); REQUIRE(nullptr != iUnit); cadet::IUnitOperation* const unit = reinterpret_cast(iUnit); diff --git a/test/GeneralRateModel2D.cpp b/test/GeneralRateModel2D.cpp index 9c863e989..58fe4c31d 100644 --- a/test/GeneralRateModel2D.cpp +++ b/test/GeneralRateModel2D.cpp @@ -239,14 +239,6 @@ TEST_CASE("GRM2D with 1 radial zone matches GRM", "[GRM],[GRM2D],[UnitOp],[Jacob REQUIRE(nullptr != mb); // Create a unit - cadet::IModel* const iUnitGrm = mb->createUnitOperation("GENERAL_RATE_MODEL", 0); - cadet::IModel* const iUnitGrm2d = mb->createUnitOperation("GENERAL_RATE_MODEL_2D", 0); - REQUIRE(nullptr != iUnitGrm); - REQUIRE(nullptr != iUnitGrm2d); - - cadet::IUnitOperation* const grm = reinterpret_cast(iUnitGrm); - cadet::IUnitOperation* const grm2d = reinterpret_cast(iUnitGrm2d); - cadet::JsonParameterProvider jpp = createColumnWithTwoCompLinearBinding("GENERAL_RATE_MODEL"); const double velocity = jpp.getDouble("VELOCITY"); const double colRadius = jpp.getDouble("COL_RADIUS"); @@ -262,6 +254,16 @@ TEST_CASE("GRM2D with 1 radial zone matches GRM", "[GRM],[GRM2D],[UnitOp],[Jacob jpp.set("LINEAR_SOLVER_BULK", "DENSE"); jpp.popScope(); + jpp.set("UNIT_TYPE", "GENERAL_RATE_MODEL"); + cadet::IModel* const iUnitGrm = mb->createUnitOperation(jpp, 0); + jpp.set("UNIT_TYPE", "GENERAL_RATE_MODEL_2D"); + cadet::IModel* const iUnitGrm2d = mb->createUnitOperation(jpp, 0); + REQUIRE(nullptr != iUnitGrm); + REQUIRE(nullptr != iUnitGrm2d); + + cadet::IUnitOperation* const grm = reinterpret_cast(iUnitGrm); + cadet::IUnitOperation* const grm2d = reinterpret_cast(iUnitGrm2d); + // Set WENO order const int wenoOrder = 3; cadet::test::column::setWenoOrder(jpp, wenoOrder); @@ -275,8 +277,8 @@ TEST_CASE("GRM2D with 1 radial zone matches GRM", "[GRM],[GRM2D],[UnitOp],[Jacob REQUIRE(grm2d->numDofs() == grm->numDofs()); - const cadet::active flowIn[] = {velocity * colPorosity * crossSectionArea}; - const cadet::active flowOut[] = {velocity * colPorosity * crossSectionArea}; + const cadet::active flowIn[] = { velocity * colPorosity * crossSectionArea }; + const cadet::active flowOut[] = { velocity * colPorosity * crossSectionArea }; grm->setFlowRates(flowIn, flowOut); grm2d->setFlowRates(flowIn, flowOut); @@ -295,13 +297,13 @@ TEST_CASE("GRM2D with 1 radial zone matches GRM", "[GRM],[GRM2D],[UnitOp],[Jacob cadet::test::util::populate(yDot.data(), [=](unsigned int idx) { return std::abs(std::sin((idx + nDof) * 0.13)) + 1e-4; }, nDof); // Setup matrices - const cadet::AdJacobianParams noAdParams{nullptr, nullptr, 0u}; - grm->notifyDiscontinuousSectionTransition(0.0, 0u, {y.data(), yDot.data()}, noAdParams); - grm2d->notifyDiscontinuousSectionTransition(0.0, 0u, {y.data(), yDot.data()}, noAdParams); + const cadet::AdJacobianParams noAdParams{ nullptr, nullptr, 0u }; + grm->notifyDiscontinuousSectionTransition(0.0, 0u, { y.data(), yDot.data() }, noAdParams); + grm2d->notifyDiscontinuousSectionTransition(0.0, 0u, { y.data(), yDot.data() }, noAdParams); // Compute Jacobian - const cadet::SimulationTime simTime{0.0, 0u}; - const cadet::ConstSimulationState simState{y.data(), yDot.data()}; + const cadet::SimulationTime simTime{ 0.0, 0u }; + const cadet::ConstSimulationState simState{ y.data(), yDot.data() }; grm->residualWithJacobian(simTime, simState, jacCol1.data(), noAdParams, tls); grm2d->residualWithJacobian(simTime, simState, jacCol2.data(), noAdParams, tls); diff --git a/test/UnitOperationTests.cpp b/test/UnitOperationTests.cpp index 4e8a66962..57d942be4 100644 --- a/test/UnitOperationTests.cpp +++ b/test/UnitOperationTests.cpp @@ -49,7 +49,7 @@ namespace unitoperation cadet::IUnitOperation* createAndConfigureUnit(const std::string& uoType, cadet::IModelBuilder& mb, cadet::JsonParameterProvider& jpp) { // Create a unit - cadet::IModel* const iUnit = mb.createUnitOperation(uoType, 0); + cadet::IModel* const iUnit = mb.createUnitOperation(jpp, 0); REQUIRE(nullptr != iUnit); cadet::IUnitOperation* const unit = reinterpret_cast(iUnit);