Fully adaptive multiresolution schemes for strongly degenerate parabolic equations with discontinuous flux
A fully adaptive finite volume multiresolution scheme for one-dimensional strongly degenerate parabolic equations with discontinuous flux is presented. The numerical scheme is based on a finite volume discretization using the approximation of Engquist-Osher for the flux and explicit time stepping. An adaptive multiresolution scheme with cell averages is then used to speed up CPU time and memory requirements. A particular feature of our scheme is the storage of the multiresolution representation of the solution in a dynamic graded tree, for sake of data compression and to facilitate navigation. Applications to traffic flow with driver reaction and a clarifier-thickener model illustrate the efficiency of this method.
Key words and phrases:Multiresolution schemes, strongly degenerate parabolic equations, thresholded wavelet transform, discontinuous flux, thresholding strategy
1.1. Scope of the paper
High resolution finite volume schemes for the approximation of discontinuous solutions to conservation laws are of at least second-order accuracy in regions where the solution is smooth and resolve discontinuities sharply and without spurious oscillations. Methods of this type include the schemes described in [1, 2, 3, 4, 5, 6]. In standard situations, the solution of a conservation law
exhibits strong variations (shocks) in small regions but behaves smoothly on the major portion of the computational domain. The multiresolution technique adaptively concentrates computational effort associated with a high resolution scheme on the regions of strong variation. It goes back to Harten  for hyperbolic equations and was used by Bihari and Harten  and Roussel et al.  for parabolic equations. Important contributions to the analysis of multiresolution methods for conservation laws include [10, 11, 12].
In this paper, we present a fully adaptive multiresolution scheme and corresponding numerical experiments for strongly degenerate parabolic equations with discontinuous flux. Specifically, we consider equations of the type
where we assume that for each , the function is piecewise smooth and Lipschitz continuous, and that is a vector of scalar parameters that are discontinuous at most at a finite number of points. On the other hand, we assume that the integrated diffusion function is Lipschitz continuous and piecewise smooth with for . We admit intervals with for all , such that (1.2) degenerates into the first-order equation
wherever . If degeneracy occurs on -intervals of positive length (and not only at isolated points), (1.2) is called strongly degenerate. Clearly, solutions of (1.2) are in general discontinuous, and need to be characterized as weak solutions along with an entropy condition. Applications of (1.2) with constant parameters include models of sedimentation-consolidation processes of particulate suspensions [13, 14], two-phase flow in porous media , and traffic flow with driver reaction [16, 17]. Applications with a discontinuous parameter vector include models of traffic flow on highways with discontinuous road surface conditions [16, 18], and a model of clarifier-thickener units used in engineering applications for the continuous solid-liquid separation of suspensions [19, 20]. In the latter application, the function models sediment compressibility; the special case , in which we fall back to (1.3), corresponds to a so-called ideal suspension of rigid spherical particles forming incompressible sediments. See  for further applications.
The novelty of the present paper is that we apply an adaptive multiresolution method to one-dimensional initial-value problems for (1.2). This equation is discretized in space by a first-order conservative finite volume scheme using the Engquist-Osher approximation, for which convergence results for our class of problems are available [16, 19, 20, 21, 22]. For time discretization an explicit Euler scheme is used. The multiresolution representation of the solution allows to introduce a locally refined mesh by thresholding of the wavelet coefficients while controlling the error of the approximation. This allows us to reduce the number of costly flux evaluations with respect to the finite volume scheme on a regular fine grid. Hence, a gain in CPU time can be obtained. Furthermore, the data are efficiently represented in a dynamic graded tree data structure, which also leads to memory compression.
1.2. Multiresolution schemes
In the following, we briefly outline the underlying ideas of multiresolution schemes for conservation laws and parabolic equations. The starting point is a conservative high-resolution finite volume discretization on a uniform mesh of (1.1) or (1.3). A multiresolution analysis of the solution with subsequent thresholding of the coefficients allows an approximation with less coefficients within a given tolerance. This allows us to reduce the number of costly flux evaluations required by the high-resolution scheme, which results in a gain of efficiency. For this purpose, either point values or cell averages of the numerical solution are defined on a hierarchical sequence of nested dyadic grids. Applying a multiresolution analysis to the solution, which can be efficiently done using the fast wavelet transform, we can construct a truncated representation by simple thresholding of the obtained coefficients. This procedure yields an efficient representation of the solution on a locally refined grid while controlling the error of the approximation.
The principle of multiresolution data representation consists in considering grid averages of the data at different resolutions from the finest to the coarsest grid, and in encoding the differences between two grids. Finally, one retains only the grid averages on the coarsest grid and the set of errors (or details) for predicting the grid averages of each resolution level in this hierarchy from those of the next coarsest one. In regions where the solution is sufficiently smooth, the multiresolution coefficients are small and can hence be neglected. Thus, data can be compressed by a thresholding or truncation operation, i.e. by setting to zero those components of the representation whose multiresolution coefficients (also called wavelet coefficients or details) are in absolute value smaller than a prescribed tolerance. Thresholding allows to control the so-called perturbation error thanks to norm equivalences. The representation of the solution in physical space corresponds to a locally refined grid.
The multiresolution analysis of the numerical solution automatically detects discontinuities, since a wavelet coefficient takes into account the regularity of a function in each position and on each scale. In [7, 23, 24] Harten explored this idea and introduced multiresolution schemes for efficiently solving hyperbolic conservations laws. Using the multiresolution representation of the solution, he devised a sensor to decide at which positions of a fine mesh the flux should be exactly evaluated, and where otherwise it can be obtained more cheaply by interpolation of pre-calculated fluxes on coarser scales. Still in the context of hyperbolic conservation laws and preserving flux evaluations for all fine grid positions, Bihari and Harten  developed a second-order adaptive switch for flux evaluations, keeping an essentially non-oscillatory (ENO) scheme where multiresolution coefficients were larger than a given tolerance, and otherwise using interpolation. In , Daubechies wavelets were used as a grid refinement strategy associated with finite difference stencils on an irregular grid for solving hyperbolic equations. Centered finite differences are used in  for approximating space derivatives on sparse point approximations (SPR) obtained by interpolating wavelet transforms. An SPR-based multiresolution WENO scheme is presented in . For parabolic PDEs a finite volume method with dynamical adaptation strategy to advance the grid was developed in .
An alternative adaption strategy could be based on local a posteriori error estimates by means of residual error computation. Results of a posteriori error estimates have been reported in the literature for elliptic problems (see ), parabolic problems (see , ) and hyperbolic problems (see ,  and the references therein), but there is not known results for strongly degenerate parabolic problems. In this sense, to compute local error estimator is not easy to be realized in practice, and we prefer to concentrate our effort in the strategy based on multiscale technique. The multiresolution strategy proposed herein for strongly degenerate parabolic equations with a discontinuous flux produces a gain in computational time and in memory. The solution is efficiently represented using a graded tree data structure and the costly fluxes are computed on the locally refined grid only. The computational efficiency of the multiresolution method is related to the data compression rate, that is, to the amount of significant information preserved after thresholding in comparison with the number of grid points of the finest mesh. Thus, efficiency is measured in terms of the compression rate and CPU time.
Finally, although we limit our treatment to one space dimension, the multiresolution scheme can be extended to higher dimensional problems in different ways. One possibility is to use higher dimensional wavelet transforms constructed by a tensor product approach, and through interpolations of the numerical divergence in the sense of cell averages from coarser to finer levels, the method of predicting values hierarchically can be extended as done in . Another possibility is to explore the splitting capability of the divergence by directions as in . Fully three-dimensional computations of flame instabilities are presented in .
1.3. Strongly degenerate parabolic equations and conservation laws with discontinuous flux
Equation (1.2) combines two independent non-standard ingredients of conservation laws: the strongly degenerate diffusion term , and the flux that depends discontinuously on the spatial position . We briefly review some recent results for equations that include either ingredient.
The basic difficulty associated with degenerate parabolic equations of the type
is that their solutions need to be defined as weak, in general discontinuous solutions along with an entropy condition to select the physically relevant weak solution. In  the existence of entropy weak solutions to an initial-boundary value problem for (1.4) in the sense of Kružkov  and Vol’pert and Hudjaev [37, 38] is shown via the vanishing viscosity method, while their uniqueness is shown by a technique due to Carrillo . The well-posedness of multi-dimensional Dirichlet initial-boundary value problems for strongly degenerate parabolic equations is shown in . Further recent contributions to the analysis of strongly degenerate parabolic equations include [41, 42, 43, 44].
Evje and Karlsen  show that explicit monotone finite difference schemes  converge to entropy solutions for the Cauchy problem for (1.4). These results are extended to several space dimensions in . The convergence of finite volume schemes for initial-boundary value problems is proved in [44, 48]. The monotone scheme used for numerical experiments in [19, 20] is the robust Engquist-Osher scheme . Thus, (1.4) admits a rigorous convergence analysis for suitable numerical schemes.
In the context of the clarifier-thickener model, the analysis of (1.2) for the case , that is, of the first-order conservation law with discontinuous flux (1.3), has been topic of a recent series of papers including [19, 50, 51], in which a rigorous mathematical (existence and uniqueness) and numerical analysis is provided. The main ingredient in these clarifier-thickener models is equation (1.3), where the (with respect to , nonconvex) flux and the discontinuous vector-valued coefficient are given functions. When is smooth, Kružkov’s theory  ensures the existence of a unique and stable entropy weak solution to (1.3). Kružkov’s theory does not apply when is discontinuous. In , a variant of Kružkov’s notion of entropy weak solution for (1.3) that accounts for the discontinuities in is introduced and existence and uniqueness (stability) of such entropy solutions in a certain functional class are proved. The existence of such solutions follows from the convergence of various numerical schemes such as front tracking , a relaxation scheme [51, 52], and upwind difference schemes .
Strongly degenerate parabolic equations with discontinuous fluxes are studied in [21, 22, 53]. In  equations like (1.2) are studied with a concave convective flux and with replaced by . Existence of an entropy weak solution is established by proving convergence of a difference scheme of the type discussed in this paper. Uniqueness and stability issues for entropy weak solutions are studied in  for a particular class of equations. These analyses are extended to the traffic and clarifier-thickener models studied herein in  and , respectively.
1.4. Time discretization, space discretization, and numerical stability
The numerical scheme for the solution of (1.2) is described in . In this work, the basic scheme is first order in time and space. We utilize a simple explicit Euler discretization in time. The spatial discretization is done by using the Engquist-Osher approximation for the convective part of the flux combined with a second-order conservative discretization of the diffusion term. For stability we need to satisfy a CFL condition requiring that in general be bounded. In some cases without diffusion (Example 2 of Section 6) we need only that be bounded.
1.5. Outline of this paper
The remainder of this paper is organized as follows. In Section 2, we briefly outline two applicative models that lead to an equation of the type (1.2), namely, a model of traffic flow with driver reaction and discontinuous road surface conditions (Section 2.1) and a clarifier-thickener model (Section 2.2). For detailed derivations of both models, we refer to  and , respectively. In Section 3, we describe the basic numerical finite volume scheme for the discretization of (1.2) on a uniform grid.
In Section 4, the conservative adaptive multiresolution discretization is introduced. Details on the numerical method and on its implementation using dynamical data structures can be found in . For the particular application to strongly degenerate parabolic equations with a flux that depends on but not on , we refer to .
The basic motivation of this approach is to accelerate a given finite volume scheme on a uniform grid without loosing accuracy. The principle of the multiresolution analysis is to represent a set of data given on a fine grid as values on a coarser grid plus a series of differences at different levels of nested dyadic grids. These differences contain the information of the solution when going from a coarse to a finer grid. An appealing feature of this data representation is that coefficients are small in regions where the solution is smooth. Applying a thresholding of small coefficients a locally refined adaptive grid is defined. The threshold is chosen in such a way to guarantee that the discretization error of the reference scheme is balanced with the accumulated thresholding error which is introduced in each time step. This yields a memory and CPU time reduction while controlling the precision of the computations. The dynamic graded tree is introduced in Section 4.1, while the multiresolution transform of a function, which is stored in the graded tree, is outlined in Section 4.2. The complete multiresolution algorithm is outlined in Section 4.3.
An error analysis, which has been adapted from Cohen et al.  and is also advanced in  for strongly degenerate parabolic equations of the type (1.4), is presented in Section 5. This error analysis motivates the choice of two parameters in the thresholding algorithm. In Section 6 we present three numerical examples, namely the traffic model (Example 1, Section 6.1), a sub-case of the clarifier-thickener model with that illustrates the application of the method to (1.3) (Example 2, Section 6.2), and the clarifier-thickener model treating a flocculated suspension, now again with a degenerate diffusion term (Example 3, Section 6.3). Numerical results, limitations and extensions of the method are discussed in Section 7.
2. Applications of strongly degenerate parabolic equations
2.1. Traffic flow with driver reaction and discontinuous road surface conditions
The classical Lighthill-Whitham-Richards (LWR) kinematic wave model [54, 55] for unidirectional traffic flow on a single-lane highway starts from the principle of “conservation of cars” for and , where is the density of cars as a function of distance and time and is the velocity of the car located at position at time . The decisive constitutive assumption of the LWR model is that is a function of only, . In other words, it is assumed that each driver instantaneously adjusts his velocity to the local car density. A common choice is , where is a maximum velocity a driver assumes on a free highway, and is a hindrance function taking into account the presence of other cars that urges each driver to adjust his speed. Thus, the flux is
The diffusively corrected kinematic wave model (DCKWM) [16, 17] extends the LWR model by a strongly degenerating diffusion term. This model incorporates a reaction time , representing drivers’ delay in their response to events, and an anticipating distance , which means that drivers adjust their velocity to the density seen an anticipating distance ahead. In fact, we adopt the equation  , where the first argument is the distance required to decelerate to full stop from speed at deceleration , and the second imposes a minimal anticipation distance, regardless of how small the velocity is. If one assumes that the effects of reaction time and anticipation are only relevant when the local car density exceeds a critical value , then the final governing equation (replacing ) of the DCKWM is the strongly degenerate parabolic equation
We assume that is chosen such that for . Consequently, the right-hand side of (2.3) vanishes on the interval , and possibly at the maximum density . Thus, the governing equation of the DCKWM model (2.3) is strongly degenerate parabolic.
Following Mochon , Bürger and Karlsen  extend the DCKWM traffic model to variable road surface conditions by replacing the coefficient in by a discontinuously varying function . However, the degenerate diffusion term models driver psychology and should therefore not depend on road surface conditions. Consequently, the new model equation for the traffic model is
For simplicity, we assume that on the major part of the highway, the maximum velocity assumes a constant value , which is also used as the value of entering the definition of in (2.4), and that there is an interval on which the maximum velocity assumes an exceptional value :
The initial-value problem for equation (2.5) with Cauchy data for is well posed , but we here insist on using a finite domain that can completely be represented by our data structure. Therefore we consider a circular road of length , the initial condition
and the periodic boundary condition
2.2. Clarifier-thickener model
The analysis of (1.4) has in part been motivated by a theory of sedimentation-consolidation processes of flocculated suspensions [13, 20], in which the unknown is the solids concentration as a function of time and depth . The particular suspension is characterized by the hindered settling function and the integrated diffusion coefficient , which models the sediment compressibility. The function is assumed to be continuous and piecewise smooth with for and for and . A typical example is
where is the settling velocity of a single particle in unbounded fluid. Moreover, we have that
Here, is the solid-fluid density difference, is the acceleration of gravity, and is the derivative of the material specific effective solid stress function . We assume that the solid particles touch each other at a critical concentration value (or gel point) , and that
This implies that for , such that also this application motivates a strongly degenerate parabolic equation (1.3). A typical function is
The extension of the one-dimensional sedimentation-consolidation equation (1.4) (if and have the interpretation given herein) to continuous sedimentation processes leads to the so-called clarifier-thickener model , see Figure 1. We consider a cylindrical vessel of constant cross-sectional area , which occupies the depth interval with and . At depth , fresh suspension of a given feed concentration is pumped into the unit at a volume rate . Within the unit, the feed flow is divided into an upwards directed and a downwards-directed bulk flow with the signed volume rates and , where conservation of suspension implies . Furthermore, we assume that the feed suspension is loaded with solids at the given feed concentration . Finally, at and , overflow and underflow pipes are provided through which the material leaves the clarifier-thickener unit. We assume that the solid and the fluid phases move at the same velocity through these pipes, so that the solid-fluid relative velocity is zero for and , which means that the term is “switched off” outside . See  for details.
We only consider vessels with a constant interior cross-sectional area and define the velocities and . Then the final clarifier-thickener model is given by (1.2) with
where we use the discontinuous parameters
We assume the initial concentration distribution
Thus, the clarifier-thickener model is specified by (1.2) with the discontinuous fluxes defined by the continuous functions given by (2.9) and (2.10), the discontinuous parameters (2.13), (2.14), and the initial condition (2.15).
3. Numerical scheme
The numerical scheme for the solution of (1.2) is essentially described in . We begin the definition of the base algorithm discretizing into cells , where with . Let , and . For we define the approximations according to
The symbols are spatial difference operators: and , and we use the Engquist-Osher flux 
Note that our pointwise discretization of , (3.2), follows the usage of [16, 20, 21, 22], but differs from that of , where is discretized by cell averages taken over the cells , where , . The important point is that in both cases, the discretization of is staggered with respect to that of the conserved quantity , and this property greatly facilitates the convergence analysis of the numerical schemes. If the discretizations were aligned (i.e., not staggered), we would have to deal with more complicated Riemann problems at cell boundaries. Further discussion of this point is provided e.g. in . Our particular choice of (3.2) (as opposed to forming cell averages) is basically its simplicity.
The space-time parameters are chosen in such way that we have the following CFL condition (see ):
which means that only must be bounded.
Let us mention that the scheme also admits a semi-implicit variant, in which the diffusion terms are evaluated at the time level . This variant has been used for numerical examples in , and its convergence for a similar equation with a convective flux that does not depend on , but which is supplemented by boundary conditions, has been proved in . The advantage of a semi-implicit scheme is that it is stable under the CFL condition (3.4), which is milder than (3.3), so that much larger time step could be used. However, a semi-implicit version involves the solution of systems of nonlinear equations for each time step, and these equations have to be solved iteratively by appropiate linearization. Since we with to keep the basic scheme as simple as possible and focus on the multiresolution device, we have decided to avoid this additional effort here. Additional complications possibly arise from the fact that we herein implement the scheme on an adaptive grid; a semni-implicit variant would, for example, generate nonlinear systems of different size in each time step. In general, implicit multiresolution schemes have been explored little so far.
4. Conservative adaptive multiresolution discretization
4.1. The graded dynamic tree
The reference standard finite volume scheme described in Section 3 yields solutions represented by vectors containing approximated cell averages on a dyadic uniform grid at time .
An important feature of our scheme is that the differences at different levels, and the solution at different levels, are always organized in a tree structure that is dynamic in the following sense: whenever an element is included in the tree, all other elements corresponding to the same spatial region in coarser resolutions, are also included. The data structure is organized as a dynamic graded tree mainly for sake of data and time compression and in particular to be able to navigate through the tree. The adaptive grid corresponds to a set of nested dyadic grids generated by refining recursively a given cell depending on the local regularity of the solution. The basis of the tree is called root. A node is an element of the tree. In one dimensional space, a parent node has two sons, and the sons of the same parent are called brothers. A given node has nearest neighbors in each direction, called nearest cousins. A node without sons is called a leaf. For the computation of the fluxes of a leaf, we need nearest cousins in each direction. If these do not exist, we create them as virtual leaves. In Figure 2 we illustrate the graded tree structure.
The nodes of the tree are the control volumes. Following , we denote by the set of indices of existing nodes, by the restriction of to the leaves, and by the restriction of to a multiresolution level , .
To estimate the cell averages of on level from those of the next finer level , we use the projection operator . This operator is exact, unique, and in our one-dimensional case is defined by
4.2. The multiresolution transform
To estimate the cell averages of a level from the ones of the immediately coarser level , we use the prediction operator . This operator gives an approximation by interpolation of at level . In contrast to the projection operator, there is an infinite number of choice for the definition of , but we impose two constraints: firstly, the prediction is local in the sense that the interpolation for a son is made from the cell averages of its parent and the s nearest cousins of its parent; and secondly, the prediction is consistent with the projection in the sense that it is conservative with respect to the coarse grid cell averages or equivalently, .
For a regular grid structure in one space dimension, we use a polynomial interpolation:
The detail is the difference between the exact and the predicted value: . Given that a parent has two sons, only one detail is independent. Then, knowledge of the cell average values of the two sons is equivalent to that of the cell average value of the father and the independent detail. Repeating this operation recursively on levels, we get the multiresolution transform on the cell average values
One of the features of this adaptive multiresolution discretization lies in the possibility to avoid considering the prediction error of the numerical flux in the update of the numerical solution, as in Harten’s original approach. This feature may be seen as an advantage in the frame of equations with discontinuous flux.
4.3. Multiresolution algorithm
Now we give a brief description of the multiresolution procedure used to solve the test problems.
Initialization of parameters: Model, FV and multiresolution parameters.
Creating the initial graded tree structure:
Create the root of the tree and compute its cell-average value.
Split the cell, compute the cell-average values in the sons and compute the corresponding details.
Apply the thresholding strategy for the splitting of the new sons.
Repeat this until all sons have details below the required tolerance .
Determine the set of leaves and virtual leaves.
Time evolution with fixed time step: Compute the discretized divergence operator for all the leaves. Performing of the space discretization is done in a locally uniform grid (regarding each leaf as a control volume of an uniform grid), so we need only those cell average values which are involved in the evaluation of the fluxes for the “edges” of the adaptive mesh formed by the leaves of the tree, i.e. we need the leaves and the nearest cousins in each direction.
Updating the tree structure:
Recalculate the values on the nodes and the virtual nodes by projection from the leaves. Compute the details in the whole tree. If the detail in a node is smaller than the prescribed tolerance, then the cell and its brothers are deletable.
If some node and all its sons are deletable, and the sons are leaves without virtual sons, then delete sons. If this node has no sons and it is not deletable and it is not at level , then create sons.
Update the values in the new sons by prediction operator from the former leaves.
END DO .
Output: Save mesh, leaves and cell-averages. Deallocate tables and plots.
With such a process we obtain high order approximation in the smooth regions and mesh refinement near discontinuities as a consequence of the polynomial exactness in the multiresolution prediction operator, even in the reference finite volume scheme is low order accurate.
For a given case of simulation, the performance of the multiresolution method can be assessed by two quantities: the data compression rate and the speed-up factor . The data compression rate is defined by
where and are numbers of points of the finest grid and of the leaves in the graded tree, respectively. Note that the data compression rate measures the memory compression at a given time of the simulation.
The speed-up factor is the ratio between the CPU time of the numerical solution obtained by the FV method and the CPU time of the numerical solution obtained by the multiresolution method:
5. Error analysis of the adaptive multiresolution scheme
The main properties of the basic finite volume scheme, i.e., its contractivity, the CFL stability condition and the order of approximation in space, allow to derive the optimal choice of the threshold parameter for the adaptive multiresolution scheme. Following the ideas introduced by Cohen et al.  and thereafter extended to parabolic equations by Roussel et al.  we decompose the global error between the cell average values of the exact solution at the level , denoted by , and those of the multiresolution computation with a maximal level , denoted by , into two errors
The first error on the right-hand side, called discretization error, is the one of the finite volume scheme on the finest grid of level . It can be bounded by
provided that is the convergence order of the finite volume scheme. The classical approach of Kuznetsov  allows to obtain for a hyperbolic scalar equation. Excepting the discontinuity due to the degeneracy, we can anticipate that the value is a pessimistic estimate of the convergence rate for our case. Unfortunately to our knowledge, no theoretical result of convergence rate for numerical schemes for strongly degenerated parabolic equations is available so far. Some numerical tests in  give , slightly over our chosen value.
For the second error, called perturbation error, Cohen et al.  assume that the details on a level are deleted when smaller than a prescribed tolerance . Under this assumption, they show that if the numerical scheme, i.e. the discrete time evolution operator is contractive in the chosen norm, and if the tolerance at the level is set to , then the difference between finite volume solution on the fine grid and the solution obtained by multiresolution accumulates in time and satisfies
where and denotes the number of time steps.
On the other hand, denoting by the size of the domain and the smallest space step, we have . Thus, according to the CFL condition (3.3), the time step must satisfy
If we want the perturbation error (5.3) to be of the same order as the discretization error (5.2), we need that . Following Cohen et al. , we define the so-called reference tolerance as . This gives
For the case (see Example 2 of Section 6), the reference tolerance must be taken as
because of the less restrictive CFL condition. To choose an acceptable value for the factor , a series of computations with different tolerances are necessary.
6. Numerical results
6.1. Example 1: Diffusively corrected kinematic traffic model with changing road surface condition.
Our numerical example for this model has been chosen in such a way that results can be compared with simulations shown in Example 5 of . The velocity function is given by (2.2) with , and , so that
We choose mph and mph. The initial density is chosen as
In Example 1, we consider an initial convoy of cars traveling on an empty road, and wish to see how the convoy passes through the reduced speed road segment. The numerical solution obtained by our method is represented in a three-dimensional plot in Figure 3 and shown at four different times in Figures 4 and 5. These figures also display the corresponding position of the leaves. For these four times, Table 1 displays the corresponding values of the speed-up factor , the compression rate , and normalized approximate errors. These errors and the speed-up factor are measured with respect to a fine grid calculation (no multiresolution) with cells. (We further comment on the behaviour of and in the discussion of Example 3.)
For this example, we take an initial dynamic graded tree, allowing multiresolution levels. We use a fixed time step determined by , thus . The prescribed tolerance is obtained from (5.4), where the constant for this example corresponds to a factor , so and the thresholding strategy is .
The errors in norm between the numerical solution obtained by our multiresolution scheme for different multiresolution levels , and the numerical solution by finite volume approximation in a uniform fine grid with control volumes, are depicted in Fig. 6. In practice, we compute the error between the numerical solution obtained by multiresolution and the projection of the the numerical solution by finite volume approximation. We also observe the same slope () between finite volume and multiresolution computation in the error of Figure 6.
6.2. Example 2: Clarifier-thickener treating an ideal suspension ().
For Example 2, we choose the same parameters as in [50, 51], so that results can be compared. In particular, we consider an ideal suspension that does not form compressible sediments, i.e., we set , so that the model considered in this example actually corresponds to the first-order equation (1.3).
We consider a clarifier-thickener unit that is initially full of water by setting . At , we start to fill up the device with feed suspension of concentration . We also consider and and we assume that the mixture leaving the unit at and is transported away at the bulk flow velocities and . The suspension is characterized by the function given by (2.9) with , and .
We use an initially graded tree with multiresolution levels and a reference tolerance of . The finest grid has control volumes and we choose a factor . Observe that the visual grid used to display Figure 7 coincides with the computational grid in direction, but in direction, only every 50th profile is plotted.
For Example 2, we use as a reference solution a fine grid computation with control volumes. Table 2 lists the behaviour of the error and the gain in computational effort and data storage for different times. Also, analogously to Example 1, we can observe in Figure 10 that the plots of the error, which is measured here for , have the same slopes.
6.3. Example 3: Clarifier-thickener treating a flocculated suspension ()
The parameter of the flux is the same as in Example 2, and the function is given by (2.12) with , and . The remaining parameters are and . Note that for (2.9) with , the function has an explicit closed-form representation, see . The reference numerical scheme is (3.1) with .
Our simulation corresponds to the choice and . The feed concentration corresponds to . Figure 11 shows the numerical solution until . In this case we consider the device with an initial concentration distribution of , and we can observe the initial stage of the fill-up process. For this example, we take an initial dynamic graded tree, allowing multiresolution levels, and for the reference tolerance we use , so .
For Example 3 we use as a reference solution a fine grid computation with control volumes. Table 3 again displays the behaviour of the error and the gain in computational effort and data storage for different times. Also, analogously to the previous examples, we observe in Figure 14 the same slope between the errors for the finite volume and multiresolution methods. This error is masured here at .
Note that in all numerical examples, the speed-up factor increases as is increased, even if the data compression rate remains constant, which approximately is the case in Table 3, or even decreases, as we see, for example, by comparing the values of for and in Table 1. The explanation of this discrepancy is that while measures the quality of performance of the multiresolution seen at the instant , the speed-up factor is referred to the total time of simulation and also includes the “overhead” required by initializing the graded tree in step (2) of the multiresolution algorithm. The initialization requires a fixed amount of CPU time, which is independent of the number of total time steps (which is porportional to , since we consider to be fixed). On the other hand, a standard FV method on a fixed grid will always require CPU time proportional to the number of time steps. This explains why even if does not change significantly, we observe an inprovement of the speed-up factor as is increased.
Before discussing our results, we comment that the standings of both applicative models are slightly different. Numerous mathematical models have been proposed for one-directional flows of vehicular traffic; reviews of this topic are given in the monographs by Helbing , Kerner  and Garavello and Piccoli , as well as in the articles by Bellomo et al. [65, 66, 67]. These and other works vividly illustrate that the number of balance equations (for the car density, velocity, and possibly other flow variables) that form a time-dependent model based on partial differential equations, as well as the algebraic structure of these equations, is a topic of current research. Fortunately, all these models are spatially one-dimensional, and a circular road with periodic boundary conditions provides a setup that is both physically meaningful (since the flow is horizontal) and easy to implement for numerical simulation. This setup, on one hand, is widely used to compare different traffic models, and on the other hand allows to assess the local influence and long-term behaviour of nonlinearities and inhomogeneities such as the ones introduced in Section 6.1.
While the traffic problem highlights the use of the scheme used herein to explore different models, the clarifier-thickener model calls for an efficient tool to perform simulations, on one hand, related to clarifier-thickener design and control [20, 68], and on the other hand, to parameter identification calculations [69, 70]. In fact, depending on the parameters, clarifier-thickener operations such as fill-up may extend over weeks and months , and require large simulation times, while the parameter identification procedures in [69, 70] proceed by solution of an adjoint problem, which needs storage of the complete solution of the previously solved direct problem. Clearly, methods that imply savings in both computational time and memory storage, such as the multiresolution scheme presented herein, are of significant practical interest for the clarifier-thickener model.
Both mathematical models considered herein exhibit three types of fronts that typically occur in solutions of (1.2), namely standard shocks (i.e., discontinuities between solution values for both of which (1.2) is hyperbolic), hyperbolic-parabolic type-change interfaces (such as the sediment level in Example 3), and stationary discontinuities located at the discontinuities of . The basic motivation for applying a finite volume multiresolution scheme is that this device is is sufficiently flexible to produce the refinement necessary to properly capture all these discontinuities, and leads to considerable gains in storage as can be seen from the sparsity of the graded trees in our numerical examples. Moreover, Figure 6 confirms that we may effectively control the perturbation error, in the sense that the error of the resulting finite volume multiresolution scheme remains of the same order as that of the finite volume scheme on a uniform grid. We recall from Section 5 that the feasibility of this control depends on an estimate of the convergence rate of the basic discretization on a uniform grid, which is an open problem for strongly degenerate parabolic equations.
Although our numerical results look promising, they still alert to some shortcomings that call for improvement. The most obvious one is the limitation of the time step according to the spatial step size of the finest grid, which can possibly be removed by using a space-time adaptive scheme such as the recent finite volume multiresolution schemes by Stiriba and Müller . On the other hand, the basic finite volume scheme accurately resolves the discontinuities of the solution sitting at the jumps of at any level of discretization; these discontinuities are not approximated by smeared transitions (as are discontinuities at positions where is smooth), see . This means that the refinement the multiresolution produces near these discontinuities, which is visible in Figures 8 and 9, and which is based on the adaptation of the refinement according to features of the solution (but not of ), is possibly unnecessary, and that a more efficient version of the present method may be feasible.
RB and MS acknowledge support by Fondecyt projects 1050728, 1070694, Fondap in Applied Mathematics, and DAAD/Conicyt Alechile project. KS and MS acknowledge support by Fondecyt of International Cooperation project 7050230. RR acknowledges support by Conicyt Fellowship.
-  Harten A (1983) High resolution schemes for hyperbolic conservation laws. J Comput Phys 49:357–393
-  Shu CW, Osher S (1989) Efficient implementation of essentially non-oscillatory shock-capturing schemes II. J Comput Phys 83:32–78
-  Shu CW (1998) Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws. In: Cockburn B, Johnson C, Shu CW, Tadmor E, Advanced Numerical Approximation of Nonlinear Hyperbolic Equations (Quarteroni A, Ed.), Lecture Notes in Mathematics vol. 1697, Springer-Verlag, Berlin (1998), 325–432
-  Kurganov A, Tadmor E (2000) New high resolution central schemes for nonlinear conservation laws and convection-diffusion equations. J Comput Phys 160:241–282
-  Nessyahu H, Tadmor E (1990) Non-oscillatory central differencing for hyperbolic conservation laws. J Comput Phys 87:408–463
-  Toro EF, Billett SJ (2000) Centered TVD schemes for hyperbolic conservation laws. IMA J Numer Anal 20:47–79
-  Harten A (1995) Multiresolution algorithms for the numerical solution of hyperbolic conservation laws. Comm Pure Appl Math 48:1305–1342
-  Bihari BL, Harten A (1995) Application of generalized wavelets: a multiresolution scheme. J Comp Appl Math 61:275–321
-  Roussel O, Schneider K, Tsigulin A, Bockhorn H (2003) A conservative fully adaptive multiresolution algorithm for parabolic PDEs. J Comput Phys 188:493–523
-  Dahmen W, Gottschlich-Müller B, Müller S (2001) Multiresolution schemes for conservation laws. Numer Math 88:399–443
-  Cohen A, Kaber SM, Müller S, Postel M (2001) Fully adaptive multiresolution finite volume schemes for conservation laws. Math Comp 72:183–225
-  Müller S (2003) Adaptive Multiscale Schemes for Conservation Laws. Lecture Notes in Computational Science and Engineering vol. 27, Springer-Verlag, Berlin.
-  Berres S, Bürger R, Karlsen KH, Tory EM (2003) Strongly degenerate parabolic-hyperbolic systems modeling polydisperse sedimentation with compression. SIAM J Appl Math 64:41–80
-  Bürger R, Evje S, Karlsen KH (2000) On strongly degenerate convection-diffusion problems modeling sedimentation-consolidation processes. J Math Anal Appl 247:517–556
-  Espedal MS, Karlsen KH (2000) Numerical solution of reservoir flow models based on large time step operator splitting methods. In: Espedal MS, Fasano A, Mikelić A (Eds.), Filtration in Porous Media and Industrial Application, Lecture Notes in Mathematics vol. 1734, Springer-Verlag, Berlin (2000), 9–77
-  Bürger R, Karlsen KH (2003) On a diffusively corrected kinematic-wave traffic model with changing road surface conditions. Math Models Meth Appl Sci 13:1767–1799
-  Nelson P (2002) Traveling-wave solutions of the diffusively corrected kinematic-wave model. Math Comp Modelling 35:561–579
-  Mochon S (1987) An analysis of the traffic on highways with changing surface conditions. Math Modelling 9:1–11
-  Bürger R, Karlsen KH, Risebro NH, Towers JD (2004) Well-posedness in and convergence of a difference scheme for continuous sedimentation in ideal clarifier-thickener units. Numer Math 97:25–65
-  Bürger R, Karlsen KH, Towers JD (2005) A mathematical model of continuous sedimentation of flocculated suspensions in clarifier-thickener units. SIAM J Appl Math 65:882–940
-  Karlsen KH, Risebro NH, Towers JD (2002) On an upwind difference scheme for degenerate parabolic convection-diffusion equations with a discontinuous coefficient. IMA J Numer Anal 22:623–664
-  Karlsen KH, Risebro NH, Towers JD (2003) stability for entropy solutions of nonlinear degenerate parabolic convection-diffusion equations with discontinuous coefficients. Skr K Nor Vid Selsk 3:1–49
-  Harten A (1994) Adaptive multiresolution schemes for shock computations. J Comput Phys 115:319–338
-  Harten A (1996) Multiresolution representation of data: a general framework. SIAM J Numer Anal 33:1205–1256
-  Jameson L (1998) A wavelet-optimized, very high order adaptive grid and order numerical method. SIAM J Sci Comput 19:1980–2013
-  Holmström M (1999) Solving hyperbolic PDEs using interpolating wavelets. SIAM J Sci Comput 21:405–420
-  Bürger R, Kozakevicius A (2007) Adaptive multiresolution WENO schemes for multi-species kinematic flow models. J Comput Phys 224:1190–1222
-  Verfürth R (1995) A review of A posteriori Estimation and Adaptative Mesh-Refinement Techniques. Advances in Numerical Mathematics. Wiley/Teubner, New York-Stuttgart.
-  Eriksson K, Hohnson C (1991) Adaptative finite element methods for parabolic problems. I. A linear model problem. SIAM J Numer Anal 28:43–77
-  Ohlberger M (2001) A posteriori error estimates for vertex centered finite volume approximations of convection-diffusion-reaction equations. M2AN Math Model Numer Anal 35:355–387
-  Tadmor E (1991) Local error estimates for discontinuous solutions of nonlinear hyperbolic equations. SIAM J Numer Anal 28:891–906
-  Cockburn B, Coquel F, LeFloch P (1994) An error estimate for finite volume methods for multidimensional conservation laws. Math Comp 63:77–103
-  Bihari BL, Harten A (1997) Multiresolution schemes for the numerical solution of 2-D conservation laws I. SIAM J Sci Comput 18:315–354
-  Chiavassa G, Donat R (2001) Point value multiscale algorithms for 2D compressive flows. SIAM J Sci Comput 23:805–823
-  Roussel O, Schneider K (2006) Numerical studies of thermodiffusive flame structures interacting with adiabatic walls using an adaptive multiresolution scheme. Combust Theory Modelling 10:273–288
-  Kružkov SN (1970) First order quasi-linear equations in several independent variables. Math USSR Sb 10:217–243
-  Vol’pert AI (1967) The spaces BV and quasilinear equations. Math USSR Sb 2:225–267
-  Vol’pert AI, Hudjaev SI (1969) Cauchy’s problem for degenerate second order quasilinear parabolic equations. Math USSR Sb 7:365–387
-  Carrillo J (1999) Entropy solutions for nonlinear degenerate problems. Arch Rat Mech Anal 147:269–361
-  Mascia C, Porretta A, Terracina A (2002) Nonhomogeneous Dirichlet problems for degenerate parabolic-hyperbolic equations. Arch Rat Mech Anal 163:87–124
-  Chen GQ, DiBenedetto E (2001) Stability of entropy solutions to the Cauchy problem for a class of nonlinear hyperbolic-parabolic equations. SIAM J Math Anal 33:751–762
-  Chen GQ, Perthame B (2003) Well-posedness for non-isotropic degenerate parabolic-hyperbolic equations. Ann Inst H Poincaré Anal Non Linéaire 20:645–668
-  Karlsen KH, Risebro NH (2003) On the uniqueness and stability of entropy solutions of nonlinear degenerate parabolic equations with rough coefficients. Discrete Contin Dyn Syst 9:1081–1104
-  Michel A, Vovelle J (2003) Entropy formulation for parabolic degenerate equations with general Dirichlet boundary conditions and application to the convergence of FV methods. SIAM J Numer Anal 41:2262–2293
-  Evje S, Karlsen KH (2000) Monotone difference approximation of solutions to degenerate convection-diffusion equations. SIAM J Numer Anal 37:1838–1860
-  Crandall MG, Majda A (1980) Monotone difference approximations for scalar conservation laws. Math Comp 34:1–21
-  Karlsen KH, Risebro NH (2001) Convergence of finite difference schemes for viscous and inviscid conservation laws with rough coefficients. MAN Math Model Numer Anal 35:239–269
-  Bürger R, Coronel A, Sepúlveda M (2006) A semi-implicit monotone difference scheme for an initial-boundary value problem of a strongly degenerate parabolic equation modelling sedimentation-consolidation processes. Math Comp 75:91–112
-  Engquist B, Osher S (1981) One-sided difference approximations for nonlinear conservation laws. Math Comp 36:321–351
-  Bürger R, Karlsen KH, Klingenberg C, Risebro NH (2003) A front tracking approach to a model of continuous sedimentation in ideal clarifier-thickener units. Nonlin Anal Real World Appl 4:457–481
-  Bürger R, Karlsen KH, Risebro NH (2005) A relaxation scheme for continuous sedimentation in ideal clarifier-thickener units. Comput Math Appl 50:993–1009
-  Karlsen KH, Klingenberg C, Risebro NH (2003) A relaxation scheme for conservation laws with a discontinuous coefficient. Math Comp 73:1235–1259
-  Karlsen KH, Risebro NH, Towers JD (2002) On a nonlinear degenerate parabolic transport-diffusion equation with a discontinuous coefficient. Electron J Diff Eqns 93:1–23
-  Lighthill MJ, Whitham GB (1955) On kinematic waves. II. A theory of traffic flow on long crowded roads. Proc Roy Soc London Ser A 229:317–345
-  Richards PI (1956) Shock waves on the highway. Oper Res 4:42–51
-  Dick AC (1966) Speed/flow relationships within an urban area. Traffic Eng Control 8:393–396
-  Greenberg H (1959) An analysis of traffic flow. Oper Res 7:79–85
-  Bürger R, Kozakevicius A, Sepúlveda M (2007) Multiresolution schemes for degenerate parabolic equations in one space dimension. Numer Meth Partial Diff Eqns 23:706–730
-  Bürger R, Ruiz R, Schneider K, Sepúlveda M (2007) Fully adaptive multiresolution schemes for strongly degenerate parabolic equations in one space dimension. Preprint 2007-03, Departamento de Ingeniería Matemática, Universidad de Concepción; submitted.
-  Kuznetsov NN (1976) Accuracy of some approximate methods for computing the weak solutions of a first order quasilinear equation. USSR Comp Math and Math Phys 16:105–119
-  Bürger R, Karlsen KH (2001) On some upwind schemes for the phenomenological sedimentation-consolidation model. J Eng Math 41:145–166
-  Helbing D (1997) Verkehrsdynamik. Springer-Verlag, Berlin.
-  Kerner BS (2004) The Physics of Traffic. Springer-Verlag, Berlin.
-  Garavello M, Piccoli, B (2006) Traffic Flow on Networks. American Institute of Mathematical Sciences, Springfield, MO, USA.
-  Bellomo N, Coscia V (2005) First order models and closure of the mass conservation equation in the mathematical theory of vehicular traffic flow. C R Mecanique 333:843–851.
-  Bellomo N, Delitalia M, Coscia V (2002) On the mathematical theory of vehicular traffic flow I. Fluid dynamic and kinetic modelling. Math Models Meth Appl Sci 12:1801–1843
-  Bellomo N, Marasco A, Romano A (2002) From the modelling of driver’s behavior to hydrodynamic models and problems of traffic flow. Nonlin Anal Real World Appl 3:339–363
-  Bürger R, Narváez A (2007) Steady-state, control, and capacity calculations for flocculated suspensions in clarifier-thickeners. Int J Mineral Process, to appear.
-  Coronel A, James F, Sepúlveda M (2003) Numerical identification of parameters for a model of sedimentation processes. Inverse Problems 19:951–972
-  Berres S, Bürger R, Coronel A, Sepúlveda M (2005) Numerical identification of parameters for a strongly degenerate convection-diffusion problem modelling centrifugation of flocculated suspensions. Appl Numer Math 52:311–337
-  Stiriba Y, Müller S (2007) Fully adaptive multiscale schemes for conservation laws employing locally varying time stepping. J Sci Comput 30:493–531