|
HyPar
1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
|
HyPar solves the following partial differential equation (PDE) using a conservative finite-difference algorithm on a Cartesian grid.
\begin{equation} \frac {\partial {\bf u}} {\partial t} = {\bf F}_{\rm hyp}\left({\bf u}\right) + {\bf F}_{\rm par}\left({\bf u}\right) + {\bf F}_{\rm sou}\left({\bf u}\right) \end{equation}
where \({\bf F}_{\rm hyp}\) is the hyperbolic term, \({\bf F}_{\rm par}\) is the parabolic term, and \({\bf F}_{\rm sou}\) is the source term. Each of these is discretized in space as described below (in the section "Spatial Discretization"), to obtain the following semi-discrete ordinary differential equation (ODE) in time:
\begin{equation} \frac {d {\bf u}} {d t} = \hat{\bf F}_{\rm hyp}\left({\bf u}\right) + \hat{\bf F}_{\rm par}\left({\bf u}\right) + \hat{\bf F}_{\rm sou}\left({\bf u}\right) \end{equation}
where \(\hat{\left(\cdot\right)}\) represents the spatially discretized terms. The governing PDE can be of any space dimension.
The hyperbolic term is of the following form:
\begin{equation} {\bf F}_{\rm hyp}\left({\bf u}\right) = -\sum_{d=0}^{D-1} \frac{\partial {\bf f}_d\left({\bf u}\right)}{\partial x_d} \end{equation}
and is discretized as:
\begin{equation} {\bf F}_{\rm hyp}\left({\bf u}\right) \approx - \sum_{d=0}^{D-1} \frac{1}{\Delta x_d} \left[ \hat{\bf f}_{d,j+1/2} - \hat{\bf f}_{d,j-1/2} \right] \end{equation}
where \(d\) is the spatial dimension, \(D\) is the total number of spatial dimensions, \(j\) denotes the grid index along \(d\). This is implemented in HyperbolicFunction().
The numerical approximation \(\hat{\bf f}_{d,j+1/2}\) of the primitive of the flux \({\bf f}_d\left({\bf u}\right)\) at the grid interfaces \(j+1/2\) is expressed as
\begin{equation} \hat{\bf f}_{d,j+1/2} = \mathcal{U}\left(\hat{\bf f}^L_{d,j+1/2},\hat{\bf f}^R_{d,j+1/2},\hat{\bf u}^L_{d,j+1/2},\hat{\bf u}^R_{d,j+1/2}\right) \end{equation}
where \(\mathcal{U}\) is an upwinding function. HyPar::Upwind points to the physical model-specific upwinding function that implements \(\mathcal{U}\), and is set by the initialization function of a specific physical model (for example Euler1DInitialize()). The physical model is specified by setting HyPar::model (read from "solver.inp" in ReadInputs()).
\(\hat{\bf f}^{L,R}_{d,j+1/2}\) are the left- and right-biased numerically interpolated values of the primitive of the flux \({\bf f}_d\left({\bf u}\right)\) at the grid interfaces and are computed using HyPar::InterpolateInterfacesHyp. They are initialized in InitializeSolvers() based on the value of HyPar::spatial_scheme_hyp (read from "solver.inp" in ReadInputs()). See interpolation.h for all the spatial discretization schemes implemented.
HyPar::HyperbolicFunction points to HyperbolicFunction().
The parabolic term can take two different forms as described below:-
No cross-derivatives: In this form, the parabolic term is of the following form:
\begin{equation} {\bf F}_{\rm par}\left({\bf u}\right) = \sum_{d=0}^{D-1} \frac {\partial^2 {\bf g}_d\left({\bf u}\right)} {\partial x_d^2} \end{equation}
where \(d\) is the spatial dimension, and \(D\) is the total number of spatial dimensions. If the parabolic function is in this form, then the physical model must specify HyPar::GFunction which must point to the function that computes \({\bf g}_d\left({\bf u}\right)\). In this case, the spatial discretization is carried out in one of two ways:-
\begin{equation} {\bf F}_{\rm par}\left({\bf u}\right) \approx \sum_{d=0}^{D-1} \frac {1}{\Delta x_d^2} \left[ \hat{\bf g}_{d,j+1/2} - \hat{\bf g}_{d,j-1/2} \right] \end{equation}
where \(j\) denotes the grid index along \(d\). \(\hat{\bf g}_{d,j+1/2}\) is the numerical approximation to the second primitive of \({\bf g}_d\left({\bf u}\right)\) and is computed using HyPar::InterpolateInterfacesPar.\begin{equation} {\bf F}_{\rm par}\left({\bf u}\right) \approx \sum_{d=0}^{D-1} \frac {1}{\Delta x_d^2} \left[ \mathcal{L}_d\left({\bf g}_d\right) \right] \end{equation}
where \(\mathcal{L}\) represents the finite-difference Laplacian operator, and is computed using HyPar::SecondDerivativePar.With cross-derivatives: In this form, the parabolic term is of the following form:
\begin{equation} {\bf F}_{\rm par}\left({\bf u}\right) = \sum_{d1=0}^{D-1}\sum_{d2=0}^{D-1} \frac {\partial^2 h_{d1,d2}\left({\bf u}\right)} {\partial x_{d1} \partial x_{d2}} \end{equation}
where \(d1,d2\) are spatial dimensions, \(D\) is the total number of spatial dimensions. If the parabolic function is in this form, then the physical model must specify HyPar::HFunction which must point to the function that computes \({\bf h}_{d1,d2}\left({\bf u}\right)\). In this case, the spatial discretization is carried out in one of two ways:-
\begin{equation} {\bf F}_{\rm par}\left({\bf u}\right) \approx \sum_{d1=0}^{D-1}\sum_{d2=0}^{D-1} \frac {1}{\Delta x_{d1} \Delta x_{d2}} \left[ \mathcal{D}_{d1}\left(\mathcal{D}_{d2}\left({\bf g}_d\right)\right) \right] \end{equation}
where \(\mathcal{D}_d\) denotes the finite-difference first derivative operator along spatial dimension \(d\), and is computed using HyPar::FirstDerivativePar.\begin{equation} {\bf F}_{\rm par}\left({\bf u}\right) \approx \sum_{d1=0}^{D-1}\sum_{d2=0,d2 \ne d1}^{D-1} \frac {1}{\Delta x_{d1} \Delta x_{d2}} \left[ \mathcal{D}_{d1}\left(\mathcal{D}_{d2}\left({\bf g}_d\right)\right) \right] + \sum_{d=0}^{D-1} \frac {1}{\Delta x_d^2} \left[ \mathcal{L}_d\left({\bf g}_d\right) \right] \end{equation}
where \(\mathcal{D}_d\) denotes the finite-difference first derivative operator along spatial dimension \(d\) (computed using HyPar::FirstDerivativePar). and \(\mathcal{L}\) represents the finite-difference Laplacian operator (computed using HyPar::SecondDerivativePar).The function pointers HyPar::InterpolateInterfacesPar, HyPar::SecondDerivativePar, and HyPar::FirstDerivativePar are set in InitializeSolvers() based on the value of HyPar::spatial_scheme_par (read from "solver.inp" in ReadInputs()). See interpolation.h, firstderivative.h, and secondderivative.h for the spatial disretization methods implemented.
Depending on which of the above forms are used, HyPar::ParabolicFunction points to either of ParabolicFunctionCons1Stage(), ParabolicFunctionNC1Stage(), ParabolicFunctionNC2Stage(), or ParabolicFunctionNC1_5Stage() (set in InitializeSolvers()).
HyPar::SourceFunction points to SourceFunction(). There is no discretization involved in general, and this function just calls the physical model-specific source function to which HyPar::SFunction points.
The semi-discrete ODE is integrated in time using explicit multi-stage time integration methods. The ODE can be written as:
\begin{equation} \frac {d {\bf u}} {d t} = {\bf F}\left({\bf u}\right) \end{equation}
where
\begin{equation} {\bf F}\left({\bf u}\right) = \hat{\bf F}_{\rm hyp}\left({\bf u}\right) + \hat{\bf F}_{\rm par}\left({\bf u}\right) + \hat{\bf F}_{\rm sou}\left({\bf u}\right) \end{equation}
The following explicit time integration methods are implemented in HyPar (see timeintegration.h):
If compiled with PETSc (https://petsc.org/release/), HyPar can use all the time integration methods and features implemented in the TS module of PETSc. See the following for relevant documentation of PETSc time integrators:
In addition to explicit time integration, the semi-discrete ODE can be solved using
Implementation: see petscinterface.h and SolvePETSc().
Implicit and IMEX time integration:
IMEX Time Integration: For implicit-explicit (IMEX) time integration, the semi-discrete ODE can be written as follows:
\begin{equation} \frac {d {\bf u}} {d t} = {\bf F}\left({\bf u}\right) + {\bf G}\left({\bf u}\right) \end{equation}
where \({\bf F}\left({\bf u}\right)\) is integrated explicitly in time and \({\bf G}\left({\bf u}\right)\) is integrated implicitly in time. The following flags (specified in the command line or in the .petscrc file) can be used to specify which of the hyperbolic, parabolic, or source terms are treated explicitly, and which are treated implicitly.
| Term | Explicit | Implicit --------— |
|---|---|---|
| Hyperbolic \(\hat{\bf F}_{\rm hyp}\left({\bf u}\right)\) | -hyperbolic_explicit | -hyperbolic_implicit |
| Parabolic \(\hat{\bf F}_{\rm par}\left({\bf u}\right)\) | -parabolic_explicit | -parabolic_implicit |
| Source \(\hat{\bf F}_{\rm sou}\left({\bf u}\right)\) | -source_explicit | -source_implicit |
If contradictory flags are specified, i.e.,
-parabolic_explicit -parabolic_implicit
the flag specifying implicit treatment takes precedence.
\begin{equation} \hat{\bf F}_{\rm hyp}\left({\bf u}\right) = \left[\hat{\bf F}_{\rm hyp}\left({\bf u}\right) - \delta\hat{\bf F}_{\rm hyp}\left({\bf u}\right)\right] + \delta\hat{\bf F}_{\rm hyp}\left({\bf u}\right) \end{equation}
the following flags can be used to specify which of these two terms are treated explicitly and which are treated implicitly.| Term | Explicit | Implicit ----— |
|---|---|---|
| \(\left[\hat{\bf F}_{\rm hyp}\left({\bf u}\right) - \delta\hat{\bf F}_{\rm hyp}\left({\bf u}\right)\right]\) | -hyperbolic_f_explicit | -hyperbolic_f_implicit |
| \(\delta\hat{\bf F}_{\rm hyp}\left({\bf u}\right)\) | -hyperbolic_df_explicit | -hyperbolic_df_implicit |